arXiv:1503.05918vl [astro-ph.CO] 19 Mar 2015 


Draft version March 23, 2015 

Preprint typeset using DTeX style emulateapj v. 5/2/11 


TOMOGRAPHY OF THE EERMALAT 7-RAY DIFFUSE EXTRAGALACTIC SIGNAL VIA CROSS-CORRELATIONS 

WITH GALAXY CATALOGS 

Jun-Qing Xia*’^, Alessandro Cuoco^’''’^, Enzo BranchinF’’'’*, and Matteo Viel®’*” 

' Key Laboratory of Particle Astrophysics, Institute of High Energy Physics, Chinese Academy of Science, P. O. Box 918-3, Beijing 100049, P. R. China. 
^ Collaborative Innovation Center of Modem Astronomy and Space Exploration, P. R. China. 

^ Department of Physics, University of Torino, via P. Giuria 1, 10125 Torino, Italy. 

^ Istituto Nazionale di Fisica Nucleare, via P. Giuria 1, 10125 Torino, Italy. 

^ Stockholm University - AlbaNova University Center, Fysikum, SE-10691, Stockholm, Sweden. 

® Dipartimento di Matematica e Fisica, Universita degli Studi “Roma Tie”, via della Vasca Navale 84,1-00146 Roma, Italy. 

^ INFN, Sezione di Roma Tre, via della Vasca Navale 84,1-00146 Roma, Italy. 

* INAF, Osservatorio Astronomico di Roma, Monte Porzio Catone, Italy. 

® INAF Osservatorio Astronomico di Trieste, Via G. B. Tiepolo 11,1-34141, Trieste, Italy, and 
INFN, Sezione di Trieste, via Valerio 2,1-34127, Trieste, Italy. 

Draft version March 23, 2015 

ABSTRACT 

Building on our previous cross-correlation analysis (Xia et al. 2011) between the isotropic 7-ray background 
(IGRB) and different tracers of the large-scale structure of the universe, we update our results using 60-months 
of data from the Large Area Telescope (LAT) on board the Fermi Gamma-ray Space Telescope (Fermi). We 
perform a cross-correlation analysis both in configuration and spherical harmonics space between the IGRB 
and objects that may trace the astrophysical sources of the IGRB; QSOs in the Sloan Digital Sky Survey-DR6, 
the SDSS-DR8 Main Galaxy Sample, Luminous Red Galaxies (LRGs) in the SDSS catalog, infrared selected 
galaxies in the Two Micron All Sky Survey (2MASS), and radio galaxies in the NRAO VLA Sky Survey 
(NVSS). The benefit of correlating the Fenni-LAJi signal with catalogs of objects at various redshifts is to 
provide tomographic information on the IGRB which is crucial to separate the various contributions and to 
clarify its origin. 

The main result is that, unlike in our previous analysis, we now observe a significant (>3.5 tr) cross¬ 
correlation signal on angular scales smaller than 1° in the NVSS, 2MASS and QSO cases and, at lower statis¬ 
tical significance (~3.0 a), with SDSS galaxies. The signal is stronger in two energy bands, E > 0.5 GeV and 
£ > 1 GeV, but also seen at £ > 10 GeV. No cross-correlation signal is detected between Fermi data and the 
LRGs. These results are robust against the choice of the statistical estimator, estimate of errors, map cleaning 
procedure and instrumental effects. 

Finally, we test the hypothesis that the IGRB observed by £ermi-LAT originates from the summed contribu¬ 
tions of three types of unresolved extragalactic sources: BL Lacertae objects (BL Lacs), Flat Spectrum Radio 
Quasars (FSRQs) and Star-Forming Galaxies (SFGs). We find that a model in which the IGRB is mainly pro¬ 
duced by SFGs (72^37 % with 2a errors), with BL Lacs and FSRQs giving a minor contribution, provides a good 
fit to the data. We also consider a possible contribution from Misaligned Active Galactic Nuclei (MAGNs), 
and we find that, depending on the details of the model and its uncertainty, they can also provide a substantial 
contribution, partly degenerate with the SFG one. 

Subject headings: cosmology: theory - cosmology: observations - cosmology: large scale structure of the 
universe - gamma rays: diffuse backgrounds 


1. INTRODUCTION 

The origin of the extragalactic 7-ray background (EGB) 
is still unknown. After its detection and early attempts 
to unveil its origin (Kraushaar et al. 1972; Fichteletal. 
1973; Mayer-Hasselwander et al. 1982; Padovani et al. 
1993; Stecker & Salamon 1996; Sreekumar et al. 1998; 
Keshet et al. 2004; Strong et al. 2004) major advances have 
recently been possible thanks to the Fermi 7-ray Space 
Telescope. Observations with the Large Area Telescope 
(£en«/-LAT) (Atwood et al. 2009) are resolving an ever 
growing number of sources, making possible to characterize 
their properties, e.g., Ajello et al. (2012, 2014), and to 
constrain their contribution to the EGB. These constraints 
are further complemented by comparing the unresolved 
EGB with semi-analytical models of different types of 
sources, e.g., Stecker & Venters (201 1); Makiya et al. (201 1); 
Dobardzic (fe Prodanovic (2014); Tamborra et al. (2014); 


Ackermann et al. (2012a); Rephaeli & Persic (2013). Thanks 
to £erm/-LAT a sizable fraction of the EGB is starting to 
be resolved (Ackermann et al. 2014a). Therefore, to avoid 
confusion, it is convenient to use a specific term for the 
unresolved part which is the quantity we want to study, 
to distinguish it from the resolved point sources either 
masked or subtracted. In the following, we indicate the unre¬ 
solved component as the Isotropic Gamma-Ray Background 
(IGRB). 

The study of the IGRB is hampered by the presence of spu¬ 
rious contributions, like the Galactic foregrounds, or bright 
point sources, that, if not properly subtracted, contaminate 
the mean signal and generate systematic errors in the anal¬ 
ysis of the true signal. These systematic uncertainties can, 
in principle, be mitigated by considering the angular corre¬ 
lation properties of the IGRB (Ackermann et al. 2012a). In 
practice, however, the auto-correlation signal is also quite 
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prone to the aforementioned systematics, which, being not 
perfectly isotropic, affect the measurements on various an¬ 
gular scales. Instead, an effective way to enhance the sig¬ 
nal and filter out systematic effects is to cross-correlate the 
IGRB with several different distributions of extragalactic ob¬ 
jects that may or may not trace the actual source of the IGRB 
but certainly do not correlate with the sources of system¬ 
atic errors. This approach has been proposed and adopted 
by, e.g., Cuoco et al. (2008), Ando & Pavlidou (2009) and 
Xiaetal. (2011). Besides the cross-correlation with cata¬ 
logs of extragalactic objects, a further possibility recently 
proposed is to cross-correlate the IGRB with weak gravita¬ 
tional lensing maps (Shirasaki et al. 2014; Fornengo & Regis 
2014; Camera et al. 2013; Fornengo et al. 2014; Camera et al. 
2014), which presents the advantage of tracing the gravita¬ 
tional potential without any bias. Furthermore, it has been re¬ 
cently shown that cross-correlation of the IGRB with galaxy 
catalogs can provide tight constraints on the dark matter anni¬ 
hilation (Ando 2014). 

Results so far have been negative since no clear correlation 
signal has been detected with a large statistical significance. 
For example, in Xia et al. (201 1) the significance of the corre¬ 
lation between Fermi data and SDSS Luminous Red Galaxies 
(LRGs) was reported to be only at the 2a confidence level. 

Nevertheless, these results have been used to set upper lim¬ 
its on the contribution of different types of potential 7 -ray 
sources such as blazars and star-forming galaxies as well as on 
the mass and cross section of WIMP dark matter candidates 
that may also contribute to the 7 -ray background through their 
self-annihilation (e.g., Ando & Komatsu (2013); Cholis et al. 
(2014)). 

The goal of this work is to extend and improve on the origi¬ 
nal study of Xia et al. (201 1) using the most recent 7 -ray maps 
obtained with the Fernu'-LAT. We estimate the two-point an¬ 
gular cross-correlation function (CCF) and the cross-angular 
power spectrum (CAPS) of the Fer»u-LAT IGRB with a vari¬ 
ety of catalogs of objects; SDSS-DR 6 quasars (Richards et al. 
2009), SDSS-DR 8 Luminous Red Galaxies (Abdalla et al. 
2008), NVSS radiogalaxies (Blake & Wall 2002) 2MASS 
galaxies (Jarrett et al. 2000) and DR 8 SDSS main sample 
galaxies (Aiharaetal. 2011). These catalogs have in com¬ 
mon: i) a large sky coverage that could allow maximizing a 
potential cross-correlation signal; ii) the fact that they have 
been already used to perform quantitative cosmological stud¬ 
ies of the Large Scale Structures (LSS). The fact that they 
contain different types of sources that span different ranges of 
redshifts is important since it increases the sensitivity of the 
cross-correlation analysis to i) the type of sources that con¬ 
tribute to the IGRB and ii) the cosmic epoch in which this 
contribution has been provided. 

In Xia et al. (201 1) we predicted that after ten years of data 
taking by the Femu'-LAT we would be able to detect a possible 
contribution to the 7 -ray background from relatively nearby 
{z < 2) sources with a confidence level of ^ 97 %. This pre¬ 
diction was quite conservative since it was based on the ex¬ 
pected Poisson noise level. Improvements i) in the model for 
the Galactic diffuse signal, ii) in the characterization of the 
instrumental point-spread function (PSF) that allows push¬ 
ing the analysis to energies lower than 1 GeV and to scales 
smaller than 1 °, and finally. Hi) the increase in the number of 
resolved sources, allows us to improve our conservative esti¬ 
mate and justifies our decision to repeat the cross-correlation 
analysis using the 5-years Fermi maps with energies as small 
as 0.5 GeV. 


Following Xia et al. (2011) we compare the results of the 
cross-correlation analysis with theoretical predictions ob¬ 
tained under the hypothesis that the diffuse 7 -ray background 
has contributions from known extragalactic sources and set 
constraints on popular candidates like galaxies with strong 
star formation activity and two types of blazars: the flat spec¬ 
trum radio quasars (FSRQ) and the BL Lacertae (BL Lac) 
objects. 

In this work we assume a flat Cold Dark Matter model 
with a cosmological constant (ACDM) with cosmologi¬ 
cal parameters = 0.0222, = 0.1189, t = 0.095, 

h = 0.678, lnl0i% = 3.097 at kp = 0.05 Mpc'^ and n, = 
0.961 that are in agreement with recent Planck results 
(Planck Collaboration et al. 2013). 

The layout of the paper is as follows: in Section 2 we briefly 
review the theoretical background of the cross-correlation 
analysis. In Section 3 we present the Fermi maps, the various 
masks and discuss the procedure adopted to eliminate the po¬ 
tential spurious contributions to the extragalactic signal. The 
maps of the angular distribution of the extragalactic objects 
that we cross-correlate with the Fermi maps are presented in 
Section 4. In Section 5 we describe the statistical estimators 
used in our cross-correlation analysis, while in Section 6 we 
test the robustness of the results to the cleaning procedure and 
to the instrument response modeling and to the data selection. 
The results are presented in Section 7, compared to model 
predictions in Section 8 , and discussed in Section 9 in which 
we also summarize our main conclusions. 


2 . THEORETICAL BACKGROUND 


Here we briefly summarize the theoretical framework 
adopted in our analysis. In this work we use the same for¬ 
malism as in Xia et al. (201 1) to which we refer the reader for 
a more thorough discussion. 

Let us consider a population of 7 -ray sources, j, with 
power-law energy spectra I{E) oc characterized by a lu¬ 
minosity function (LE) ^j(Lj,T,E,z) in which we highlight 
the explicit dependence on the observed 7 -ray energy E, the 
rest-frame luminosity of the sources (generally expressed 
in erg s“'), cosmological redshift z, and photon index r;> 1 . 
The contribution of this population to the differential energy 
flux is; 


dE 



pLMAxiz) pTm 
/LmIN 


<l>y(L-j,,r,(l +z)E,z)L^dL^dT 


dz 

{\+z)H(z) 

( 1 ) 


where H{z) = /Fo[(l -I-z)^Hm + Ha] represents the expansion 
history in the assumed cosmological model and ( 1 -l-z) ac¬ 
counts for the cosmological redshift. All sources along the 
line of sight contribute to the integral over z- 
The integration over is performed within a finite lumi¬ 
nosity range. We set the upper value equal to 


Tmax(z) = AT : dl ( z ) Sv,M (2) 


where di is the luminosity distance in the adopted cosmol¬ 
ogy and 5iini is the (energy) flux detection limit. In general, 
the flux detection threshold depends on the power-law in¬ 
dex. This dependence is strong for the photon flux and much 
weaker for the energy flux. Eor this reason, in this work we 
shall ignore the correlation between and Tj. This im¬ 
plies that resolved sources are excluded and that the integral 
(Eq. 1) has contributions only from unresolved sources. The 
lower integration limit, Lmin, is taken from recent literature 
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TABLE 1 

Parameters of the LDDE LEs taken from Ajello et al. 2012 for FSRQs and Ajello et al. 2014 for BE Lacs. 


Model 

A“ 

71 


72 

z* 

Pi 

T 

P2 

a 


d 

<J 

BLLacsl LDDE 

3.39 X lO'* 

0.27 

0.28 

1.86 

1.34 

2.24 

4.92 

-7.37 

4.53 X 10^2 

2.10 

6.46 X 10-2 

0.26 

BLLacs2 LDDE 

9.20 X 10- 

1.12 

2.43 

3.71 

1.67 

4.50 

0.0 

-12.88 

4.46 X 10-2 

2.12 

6.04 X 10-2 

0.26 

FSRQ LDDE 

3.06 X lb* 

0.21 

0.84 

1.58 

1.47 

7.35 

0.0 

-6.5! 

0.21 

2.44 

0.0 

0.18 


a In units of 10 Mpc ^ erg * s. b In units of 10^* erg s 


for specific source classes as we shall discuss in the next sec¬ 
tion. We note that the integral converges and setting Lmin 
to much smaller values has very little effect on the final re¬ 
sults, i.e., our results are robust against the value of Lmin- The 
choice of S\im depends on the 7 -ray source catalog used to 
mask resolved point sources. In the following we will use the 
2FGL (Nolan et al. 2012) source catalog and a preliminary 
version of the 3FGL (Acero et al. 2015) catalog. Typical val¬ 
ues of the source detection thresholds (in units of integrated 
energy flux above 100 MeV) are 5 x 10“'^ erg cm“^s“' and 
2.5 X 10“'^ erg cm“^s“' respectively for the 2FGL and 3FGL 
catalog (Acero et al. 2015), with the lower threshold of the 
3FGL catalog which is in part due to the larger dataset used 
(4 years vs. 2 years) and in part to the improved character¬ 
ization of the response of the LAT. In practice, however, the 
luminosity density p^(z) that we shall use to characterize the 
contribution of a given type of sources to the energy flux is 
weakly dependent on the detection threshold, i.e., similar re¬ 
sults are obtained with the 3FGL and 2FGL thresholds, even 
if the former has a deeper reach in flux than the latter. This 
reflects the fact that below 100 GeV the bulk of the EGB is 
still unresolved in both the 2FGL and 3FGL catalogs. For this 
same reason, the energy density is insensitive to the precise 
modeling of the detection efficiency, which is not exactly a 
step function but represents rather a smooth transition in flux 
between zero and full efficiency. 

In the integration over F we assume that the intrinsic distri¬ 
bution of photon indices is a Gaussian, which implies that for 
a given redshift z and luminosity the LF has the form 

$(L-y,z,r) oc e 2-2 , ( 3 ) 

where p and a are, respectively, the mean and dispersion of 
the distribution. The mean is allowed to be a function of the 
source luminosity (expressed in units of 10 "^*^ erg s“'): 

p(L^) = /r* + /3 X (logio(L^)-46). (4) 

Since in the luminosity range of interest cr <C /i, as we will 
see, then we can approximate the F distribution with a Dirac 
delta centred on F = p. With this approximation the integrated 
flux Ij{> E) can be expressed as 

/•“ dh cL^-r^ f 


thresholds: /(> E = 0.5 GeV), /(>£ = ! GeV) and /(> E = 
10 GeV). 

Variations in the number density of unresolved sources, 
n^(z,x), are responsible for the local fluctuation in the 7 -ray 
luminosity density, pj(z,x) and, therefore, in the integrated 
7 -ray flux. If the luminosity is proportional to the number 
of sources then the two fluctuations in and p^ are related 
through 




P 7 W 


n^z) 


where the mean number density of sources is = f ^(L)dE 
We further assume a linear mapping, called biasing, between 
mass density p„, and the number density of objects: 


S„^(z,x) = b^(z)S,„(z,x) = (8) 

Pm(z) 

where we allow for a redshift-dependent bias parameter, 
b^iz). 

Fluctuations in the integrated flux along the generic direc¬ 
tion n can be obtained from (5) and the linear biasing pre¬ 
scription ( 8 ): 


Slin) = 


/(>£’,n)-/(>£) _ / p^iz)byiz)6i„(z,x)dz 
/(>£) “ /(>£) 


(9) 


where I = /(> E) indicates the 7 -ray mean flux and /(n) = 
/(> Ljn) is the energy flux along the generic direction n. 

Our goal is to investigate the cross-correlation between the 
diffuse IGRB maps and the sky-projected spatial distribu¬ 
tion of different types of extragalactic objects. The angular 
cross-power spectrum between the extragalactic background 
/(> £,n) and the fluctuation of discrete sources, j, can be ex¬ 
pressed as 


= ^Jk^P(k)[G‘i(k)][Gi(k)]dk , (10) 

where P(k) is the power spectrum of density fluctuations, I is 
the multipole of the spherical harmonics expansion and the 
functions G{k) specify the contribution of each held to the 
cross-correlation signal. More specifically, the contribution 
of the IGRB is given by 


where 

ri.MAx(z) ..y/ijiL-i) 

p^(z)= / $,(L^,z)L^^^-- dL^ (6) 

4l„in 

is the the mean luminosity density at z and $j(L-y,z) = 
4>/(L-y,z,r = p{Lj)). In this paper we deal with maps of pho¬ 
ton counts rather than energy flux; the photon flux (above 
energy E) being simply (2-rj)/(l - Fy) x Ij{> E)/E. We 
will consider maps of integrated flux above three energy 


G\{k) = J p^(z)b^(z)D(z)j,[kx(z)]dz (11) 

where ji [kx(z)] are spherical Bessel functions, D{z) is the lin¬ 
ear growth factor of density fluctuations and x(z) is the co¬ 
moving distance to redshift z. The analogous quantity for the 
number density fluctuations in a population of discrete objects 
is 

Gi(k) = J ^^bj(z)D(z)j,[kx(z)]dz, (12) 
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(e.g., Persic et al. 1989). In our case, we are required to 
model; i) the correlation properties of the underlying mass 
density field; ii) its relation with discrete tracers, i.e., the 
biasing prescription; Hi) the mean IGRB flux. To model 
the cross-correlation signal we consider the cosmologically 
evolving mass density power spectrum, Pik,z), obtained from 
the public code CAMB (Lewis et al. 2000) for the linear part 
and the Halofit built-in routine for non-linear evolution 
(Smith et al. 2003). In addition, we use the linear growth fac¬ 
tor D(z) and the comoving distance x(z) appropriate for the 
cosmological model adopted. 

To model the bias and the mean IGRB signal we need to 
select a class of objects which are likely to contribute to the 
IGRB and specify the energy spectrum, luminosity function 
and fraction of IGRB contributed by the source, fj. Here we 
consider three different candidates: FSRQs, BL Lacs and star¬ 
forming galaxies (SFGs) 


Fig. 1.— Integrated 7 -ray flux per logarithmic redshift bin dl(> E)/dlnz 
as a function of z for three different source classes: FSRQs (red, dashed), 
BL Lacs (black or magenta, continuous) and star-forming galaxies (blue or 
green, dot-dashed). 

where dN{z)/dz is the redshift distribution of the objects. 
Here we make the hypothesis that these objects, which do not 
necessarily coincide with the sources of the IGRB, also trace 
the underlying mass density held modulo some z-dependent 
linear bias parameter hj(z). 

Note that in the cross-correlation we consider the integrated 
flux /(> E,n) rather than its dimensionless analogous 6l{n) 
given in Eq. 9. With this choice the cross-correlation sig¬ 
nal is robust to any spurious monopole term arising from 
an incorrect subtraction of the model Galactic diffuse signal 
or charged particle contamination. One implication of this 
choice is that our model cross-correlation signal ( 10 ) is de¬ 
pendent on the mean integrated flux, /(> E). For this reason, 
and to account for uncertainties in the estimate of the mean 
IGRB signal, we allow for some freedom in the normalization 
of the luminosity function of the putative 7 -ray sources and, 
accordingly, add an additional free parameter in the model. 

In this work we also estimate the angular two-point cross¬ 
correlation function of the flux maps and discrete object cat¬ 
alogs which is simply the Legendre transform of the angular 
power spectrum 

(/(> E,ni)Sj(n2)) = ^^C‘'^P,[cos(0)] , (13) 

/ 

where Pi[x] are the Legendre polynomials and 9 is the sepa¬ 
ration angle between directions ni and n 2 . The angular two- 
point correlation function and power spectrum are two ways 
of expressing the same information. However, in practice, the 
two statistics are somewhat complementary as they probe dif¬ 
ferent scales with different efficiency and their respective es¬ 
timators are prone to different types of biases. For this reason 
we shall compute both quantities. 

2.1. Modeling the mean flux and the cross-correlation 
signal 

One of the aims of this work is to compare the measured 
cross-correlation signal with model predictions obtained un¬ 
der the assumption that some specific type of unresolved 
sources contributes to the IGRB. We note that even auto¬ 
correlation studies can provide constraints on the nature and 
spatial clustering of the sources contributing to the signal 


1. FSRQs are a type of AGN (blazars) with a relativis¬ 
tic jet pointing close to the line of sight. Ajello et al. 
( 2012 ) have recently determined the 7 -ray luminosity 
function of these objects which they have parametrized 
in the framework of a Luminosity Dependent Density 
Evolution (LDDE) model; 


4>(L^,z = 0 ,r) 


A 

In(lO)L-), 



(14) 

The term in parentheses, a smoothly-joined double 
power-law function, represents the luminosity function 
of the local FSRQs and the exponential term is the 
same photon index distribution as Eq. 3. In the LDDE 
model the luminosity function at the redshift z can be 
expressed as 


$(L.^,z,r) = ^(Lj,z = 0,r) X e{z,Lfl ), 


where 
e{z,Lj) = 
with 


1+z / i+z 

l+Zc(L-,)) ^Vl+2c(T7)y 


pi(L^) = p*-\-TX (logio(L-y)-46) 

and 

z,(L^) = z:-(L-,/104V. 


(15) 


(16) 


(17) 

(18) 


Zc represents the luminosity-dependent redshift at 
which the evolution changes from positive to nega¬ 
tive and z* is the evolutionary peak for an object with 
a luminosity of 10^^ erg s“'. This LDDE luminosity 
function model is specified by the 12 parameters listed 
in Table 1 with the particular values determined by 
Ajello et al. (2012) by fitting 7 -ray data. In the fit, 
the authors have set /? = t = 0 , i.e., they have assumed 
that neither the overall shape of the luminosity function 
nor the spectral index depend on the luminosity of the 
sources, Lj. Note that the evolutionary term e(z,Lry) 
in Eq. 16 is not equal to unity at z = 0. To derive the 
density p.y(z) in Eq. 6 required to calculate the corre¬ 
lations, we set Lmin = 10 '*"^ erg s“' as recommended 
in Ajello et al. (2012), although, as already explained. 
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choosing a lower value or even zero does not signifi¬ 
cantly affect the results. 

The parameters of the luminosity function uniquely de¬ 
termine the contribution of FSRQs to the mean diffuse 
IGRB signal. However, as anticipated, in this work we 
prefer to keep the normalisation of the luminosity func¬ 
tion free. This additional degree of freedom is meant to 
absorb experimental errors in the measurement of the 
mean diffuse IGRB signal and uncertainties in model¬ 
ing the clustering of the sources. The resulting redshift 
distribution of the 7 -ray flux contributed by FSRQs 
shown in Fig. 1 is rather broad and peaks at z ^ 0.5. 

The last ingredient of the model is the bias of the 
sources. Here we adopt the redshift-dependent AGN 
bias proposed by Bonoli et al. (2009) in the framework 
of the semi-analytic models of AGN-black holes co¬ 
evolution: Z7fSRQ(z) = 0.42-1-0.04(1 -I- z)-i-0.25(1 -I-z)^. To 
test the robustness of our results to the biasing scheme 
we have considered two alternative bias models: i) the 
rather unphysical case of a constant bias bfSRQ = 1.04 
obtained by considering it equal to the bias model of 
Bonoli et al. (2009) estimated at z = 0.5. ii) a linear, z- 
dependent model in which the bias is set equal to that of 
a 10 '^ M 0 halo. This latter choice reproduces the bias 
of X-ray selected AGN estimated by Koutoulidis et al. 
(2013) and represents an upper limit, since the bias 
of optically selected AGN is matched by the bias of 
IO'^Mq halos. 

2. BL Lacs. These sources represent a different sub-class 
of blazars with, on average, lower luminosities than FS¬ 
RQs. To model their luminosity function we adopt the 
LDDE functional form as in Eq. 14 and set the free pa¬ 
rameters according to the best fit to 7 -ray data obtained 
by Ajello et al. (2014). We consider two possibilities 
corresponding to the two sets of parameters in Table 1 . 
In the first case, that we dub BLLacsl, the authors let 
all parameters free to vary. In the model BLLacs2, they 
instead set t = 0 , i.e., basically switching off the de¬ 
pendence on the luminosity of the pi index in the evo¬ 
lutionary term. It should be stressed that the BLLacsl 
model represents a better fit to the data in Ajello et al. 
(2014). Nonetheless, we consider also model BLLacs2 
to study the robustness of the interpretation of the cross¬ 
correlation in terms of BL Lacs. In accordance with 
Ajello et al. (2014) we set Lmin = 7 x 10"^^ erg s“' for 
the calculation of the p-y(z) integral. 

The redshift distribution of the 7 -ray flux of BL Lacs 
for the two models considered is shown in Lig. 1. In 
both cases the distribution is rather broad. However, 
in BLLacsl the distribution is more local and peaks 
at z ^ 0.1, whereas in BLLacs2 the peak is at much 
higher redshift (z ~ 0.7) although with a significant tail 
at low z. We assume that the biasing of BL Lacs is 
the same in both models and equal to that of LSRQs, 
i.e., bBLLac = bfSRQ. The robustness of the results on 
this choice is also tested using the same alternative bias 
models considered for LSRQs. 

3. SFGs. As reference and comparison with our previous 
work (Xia et al. 201 1) we adopt the phenomenological 
model of Lields et al. (2010) in which the 7 -ray emis¬ 
sion in SLGs over cosmic time is proportional to their 



Fig. 2.— Comparison between the 7 -ray flux per logarithmic redshift bin 
dli> E)/dlnz. 

Star-formation rate and the gas mass-fraction, both nor¬ 
malised to the present values in the Milky Way (MW). 
The energy spectrum is assumed to be similar to the 
one observed in the MW which can be modeled ap¬ 
proximately as a broken power law with photon indexes 
F ~ 1.5 and F ~ 2.475 above and below ^ 500 MeV, 
respectively. The contribution to the IGRB, shown in 
Fig. 1, is spread over a wide redshift range and peaks at 
z ~ 1 , with a high-redshift tail more extended than that 
of BL Lacs and LSRQs. We dub this model SLGsl. 

Lor the sake of completeness we consider also 
a second model (SLGs2), originally proposed by 
Ackermann et al. (2012a) and recently revised by 
Tamborraet al. (2014). The ingredients of this model 
are the SLG LL obtained from infrared observations 
(Rodighiero et al. 2010; Gruppioni et al. 2013) and the 
empirical relation between luminosity in the infrared 
and 7 -ray bands calibrated using a samples of lo¬ 
cal galaxies observed in both bands and assumed 
to be valid at all redshifts (Ackermann et al. 2012a). 
The most recent IR observations of Gruppioni et al. 
(2013) have enabled the measurement of the LLs of 
different sub-populations of galaxies. Specifically, 
Gruppioni et al. (2013) subdivide the infrared galax¬ 
ies into normal spiral galaxies (SP), starburst galaxies 
(SB), and galaxies hosting an AGN but whose infrared 
emission is still dominated by star-forming activity (SL- 
AGN), and provide the LLs separately for each sub¬ 
population. We model the 7 -ray emission separately 
for the three populations by assuming a power-law en¬ 
ergy spectrum with F = 2.475 for SP and SF-AGN and 
of r = 2.2 for SB (Ackermann et al. 2012a). To test if 
the redshift distribution of the 7 -ray emission is robust 
with respect to the assumed spectral shape of galaxies 
we use the LF of the whole population of infrared ob¬ 
jects from Gruppioni et al. (2013). Assuming a single 
global index F = 2.475 for this population does not sig¬ 
nificantly modify the results. It has to be noted that 
besides the contribution from SP, SB and SF-AGN this 
global LF contains also a contribution from pure AGNs, 
which is, however, very subdominant. The IGRB con¬ 
tribution of SP H- SB H- SF-AGN as a function of redshift 
is plotted in Fig. 1. The distribution is significantly dif- 
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ferent than in model SFGsl since in this case the IGRB 
is mostly contributed by low-redshift galaxies. The dis¬ 
crepancy between the two distributions reflects a funda¬ 
mental difference between the two models. Both mod¬ 
els use the luminosity density in the infrared band, in 
Fields et al. (2010) assuming that is a good tracer of the 
star formation history (which is a common assumption, 
see, for example, the discussion in Rodighiero et al. 
(2010)), and in Ackermann et al. (2012a) assuming it is 
a tracer of the 7 -ray LF itself. Flowever, model SFGsl 
further contemplates the possibility of a time-dependent 
gas quenching that reduces the 7 -ray emission at low 
redshift. As a result we expect that the two models 
predict very different cross-correlation signals, despite 
having similar 7 -ray luminosities integrated over red- 
shift (Tamborra et al. 2014). 

Finally, based on the observations in Afshordi et al. 
(2004) and theoretical arguments (Wilman et al. 2008) 
we assume for both models that SFGs trace the under¬ 
lying mass density held with no bias (i.e., bspc = !)• 
We also consider an alternative case in which the bias 
of the SFGs is set equal to that of a Milky Way-sized 
haloof IO'^Mq. 

2.2. MisalignedAGNs 

Another type of source that potentially contributes to the 
IGRB is misaligned AGNs (MAGN) (Di Mauroet al. 2014; 
Inoue 201 1). Similarly to the case of SFGs, too few MAGNs 
have been detected in 7 rays to determine their LF directly in 
this band. Instead, the 7 -ray LF is inferred from that mea¬ 
sured in some other band by exploiting the observed relation 
between the luminosities in the two bands. For MAGN this 
is done by considering the radio band, i.e., by using the ob¬ 
served MAGN radio LF and the radio to 7 luminosity relation. 
The latter relation is calibrated on a local sample of objects 
for which observations in both bands are available and then 
extrapolated at all redshifts. 

We did follow this procedure and used both the radio lu¬ 
minosity function of Willott et al. (2001) and the radio 7 -ray 
luminosity relation derived in Di Mauro et al. (2014) to ob¬ 
tain the MAGN contribution to the IGRB at all redshifts. The 
result is shown in Fig. 2, which is the analogous to Fig. 1 but 
featuring only the MAGN and the SFGsl models. The main 
point here is the similarity between the two distributions that 
peak at z ~ 2 and extend to much higher redshift. We also note 
that while the SFG distribution is smooth, the MAGN exhibits 
a feature at z ~ 1. This reflects the composite nature of the 
radio-band LF that receives contributions from different types 
of objects with sharp breaks in their redshift distributions. 

For the scope of our analysis, the fact that the two redshift 
distributions are similar implies a potential degeneracy in the 
model predictions. This degeneracy can be broken only if the 
linear bias parameters of the two populations, upon which the 
amplitude of the clustering depends, are different. In fact, we 
do expect that the bias of the MAGN, which are typically as¬ 
sociated with massive dark matter halos, is higher than that of 
the SFGsl and, as a consequence, the MAGN bias is higher 
than that of SFGs at all redshifts (see e.g., Wilman et al. 
(2008), Lindsay et al. (2014) and references therein). This to¬ 
gether with the fact that the bias is an increasing function of 
redshift and that both populations contribute to the IGRB out 
to high redshifts make it possible, in principle, to discrimi¬ 
nate between the two populations through a cross-correlation 


analysis. In practice, however, the bias of both populations is 
ill-constrained by present observations. 

For this reason, in the present analysis, we focus on 
those objects for which the contribution to the IGRB and its 
anisotropy are more robust to the uncertainties in the bias. In 
fact, we have considered blazars, since their dl{> £’)/^f Inz is 
suppressed at high redshift where the bias of the objects is ex¬ 
pected to increase, and SFG, since their bias is expected to be 
close to unity even at high redshifts. For the very same reason 
we have decided not to include MAGN in our model; they 
can be found at very high redshift and their bias is large and 
rapidly increases with redshift. Of course this does not mean 
that MAGNs do not contribute to the correlation signal. Only 
that our analysis will not be able to discriminate between SFG 
and MAGN contributions. Therefore it should be kept in mind 
that the SFG contribution, which we will study in the follow¬ 
ing, may actually include a MAGN signal as well and that our 
model cross-correlation signal, entirely based on blazars and 
SFGs, could be underestimated at high redshifts. 

2.3. Further theoretical contributions 

All of our models assume that the discrete sources sample, 
according to a deterministic bias relation, the underlying mass 
density field whose two-point clustering properties can be 
quantified by a nonlinear power spectrum that can be obtained 
using CAMB + Halofit. This is an approximation that ig¬ 
nores the presence of substructures within virialized halos 
and, consequently, underestimates the power on small scales. 
In the language of the halo model (Cooray & Sheth 2002), we 
underestimate the 1-halo term contribution to the correlation 
signal. Theoretical modeling of this term is challenging since 
it depends on several characteristics both of the catalog and 
of the 7 -ray emitting population which are quite uncertain. 
For example, it is necessary to specify how the galaxies of a 
catalog are distributed on average within a DM halo of given 
mass, and, analogously, how 7 -ray sources of different lumi¬ 
nosities populate the same DM halo. Both quantities can be 
modeled, but within large uncertainties. Further details are 
discussed, for example, in Ando (2014), Camera et al. (2014) 
and Ando et al. (2007). On the contrary, the shape of the 1- 
halo term is easier to model. Assuming point-like DM halos, 
this term would simply be a constant in multipole space, or a 
delta-like term at 0 = 0° in the CCF Small distortions from a 
constant arise considering the finite extent of DM halos, and 
are typically important at very high multipoles (> 1000 ) for 
halos at low redshift (< 0.1). In the following, we will thus 
adopt a phenomenological approach and model this term as 
constant in (. with a free normalization, and check against the 
data if there is a preference for the inclusion of a 1-halo cor¬ 
relation. 

A second assumption of our model is that the sources re¬ 
sponsible for the 7 -ray signal and the various astrophysical 
sources that we cross-correlate trace the same mass density 
with different bias relations, i.e., they are not the same pop¬ 
ulation. There is, however, the possibility that they are the 
same population seen at two different wavelengths (or two 
populations that largely overlap each other). In this case one 
would expect a strong cross-correlation peak in the CCF at 
0 = 0 ° corresponding, again, to a constant in multipole space. 
With enough angular resolution, this possibility could be dis¬ 
tinguished, in principle, from a pure 1-halo term, due to the 
distortion induced by the finite extent of the DM halo in the 
latter (with the possible caveat of a positive correlation signal 
at 0 > 0 ° that may arise even in this case when the emission 
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in two different bands originates from two separate regions 
of the same object, like, for example, possibly in the case of 
the nucleus and the jet of an AGN). Typically, however, the 
LAT PSF is too large to allow discriminating between the two 
cases. We will thus model both terms as constant in multipole 
space, and consider their sum as a single contribution whose 
presence will be tested in the data. We will indicate these 
contributions collectively as l-halo-\\]s£ terms. 

3. F£/?M/-LATMAPS 

In this section we describe the IGRB maps obtained from 5 
years of Fermi-LAT data taking and the masks and templates 
used to subtract contributions from i) 7 -ray resolved sources, 
ii) Galactic diffuse emission due to interactions of cosmic rays 
with the interstellar medium and in) additional Galactic emis¬ 
sion located at high Galactic latitude in prominent structures 
such as the Fermi Bubbles (Su et al. 2010; Ackermann et al. 
2014c) and Loop I (Casandjian et al. 2009). The validity of 
the masking procedure, its effectiveness in removing fore¬ 
ground and resolved source contributions to the IGRB signal, 
and its impact on cross-correlation analysis are assessed in 
Section 6 . 

Ter»u-LAT is the primary instrument onboard the Fermi 
Gamma-ray Space Telescope launched in June 2008 
(Atwood et al. 2009). It is a pair-conversion telescope cover¬ 
ing the energy range between 20 MeV and > 300 GeV. Due to 
its excellent angular resolution 0.1° above 10 GeV), very 
large field of view (~ 2.4 sr), and very efficient rejection of 
background from charged particles, it is the best experiment 
to investigate the nature of the IGRB in the GeV energy range. 

For our analysis we have used 60 months of data from Au¬ 
gust 4th 2008 to August 4th 2013. More specifically, we 
have used the P7REP_CLEAN event selection' in order to 
ensure a low level of cosmic-ray (CR) background contam¬ 
ination. Further, to greatly reduce the contamination from 
the bright Earth limb emission we exclude photons detected 
i) with measured zenith angle larger than 100 ° or ii) when 
the rocking angle of the LAT with respect to the zenith was 
larger than 52°. In order to generate the final flux maps we 
have produced the corresponding exposure maps using the 
standard routines from the LAT Science Tools^ version 09- 
32-05, using the latest recommended P7REP_CLEAN_V15 
instrument response functions (IRFs). We use both back- 
converting and front-converting events and we checked the 
robustness of the results using either data subsample (see Sec¬ 
tion 6.3). The GaRDiAn package (Ackermann et al. 2012b, 
2009) was then used to pixelize both photon count and ex¬ 
posure maps in HEALPix^ format (Gorski et al. 2005). The 
maps contain Vpix = 3 145728 pixels with an angular size of 
~0.11°x0.11° corresponding to the HEALPix resolution pa¬ 
rameter Wide = 512. Finally, the flux maps are obtained by 
dividing the count maps by exposure maps in three energy 
ranges; E > 500 MeV, E > 1 GeV and £ > 10 GeV 

To reduce the impact of the Galactic emission on our anal¬ 
ysis focused on the IGRB, we apply a Galactic latitude cut 
1^1 > 30° in order to mask the bright emission along the 
Galactic plane. Moreover, we also exclude the region as¬ 
sociated to the Fermi Bubbles and the Loop I structure. In 

’ See http://www.slac.stanford.edu/exp/glast/groups/canda/lat_Performance.htm 
for a definition of the P7REP and P7 event selections and their characteris¬ 
tics. 

^ http://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/ 

^ http://healpix.jpl.nasa.gov/ 


Fermi-2MASS CCF E>1 GeV 



Fig . 3.— Comparison of Fermi-2MASS E > 1 GeV CCF for two different 
Galactic foreground models. 

Xia et al. (201 1) we have experimented with different latitude 
cuts and found that \b\> 30° represents the best compromise 
between pixels statistics and Galactic contamination. The cor¬ 
responding mask, obtained from the tabulated contours of the 
Fermi Bubbles given in Su et al. (2010) is shown in the bot¬ 
tom panel of Fig. 4 as the bulge-like central region together 
with the horizontal strip mask corresponding to the latitude 
cut. The mask also features a number of smaller holes placed 
at the position of all resolved sources in a preliminary version 
of the 3FGL catalog. In the £ > 1 GeV maps all sources are 
masked out with a disk of 1° radius. For £ > 0.5 GeV we 
used larger disks of 2° but only for the 500 sources with the 
highest integrated flux above 100 MeV in the catalog, while 
the remaining ones are still masked with disks of 1°. Fi¬ 
nally, to exclude the contribution from the Small and Large 
Magellanic Clouds, which are more extended, we have used 
two larger circles with 3° radius. To test the robustness of 
our results on the subtraction of resolved sources we have 
also built a similar mask using the previous 2FGL catalog 
(gll_psc_v08 . fit^ ). Further details are given in Sec¬ 
tion 6.2. When cross-correlating with a given galaxy catalog, 
the mask specific to that catalog is further employed. The 
masks for the catalogs we use can be seen in Xia et al. (201 1). 

Although we select a part of the sky at high Galactic lati¬ 
tude, the Galactic diffuse emission in this region is still sig¬ 
nificant and needs to be removed. For this purpose, and to 
check the robustness to this correction, we use two models of 
Galactic diffuse emission: ring_2year_P76_v0 . fits^ 
and gll_iem_v05_revl. f it^, which we subtract from 
the observed emission to obtain the cleaned maps. Both mod¬ 
els are based on a fit of the LAT data to templates of the HI 
and CO gas distribution in the Galaxy as well as on an inverse 
Compton model obtained with the GALPROP code®. The first 
model ring_2Year_P7 6_vO . fits is tuned to 2 years of 
P 7 data and further uses uniform flat patches to model some 
features of the diffuse sky such as the Fermi Bubbles and 
Loop I. The model gll_iem_v05_revl. fit is based on 
4 years of P7REP data and adopts an alternative procedure to 
account for residual diffuse emission involving templates of 

^ http://fermi.gsfc.nasa.gov/ssc/data/access/lat/2yr_catalog/ 

^ http://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html 

® A more detailed description can be found at 
http://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html 
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residual emission obtained in early stages of model fitting to 
construct the final model. Since part of the data has been in¬ 
troduced in some form as an a posteriori model component^ 
(Ackermannet al. 2014b), the gll_iem_v05_revl. fit 
model is not recommended for diffuse emission studies. How¬ 
ever, since these additional templates only affects regions that 
are within our masked areas, we can safely use also this model 
to test the robustness of the results with respect to the mod¬ 
eling of the Galactic emission. Indeed, the comparison be¬ 
tween the two models shows actually that they are very sim¬ 
ilar in our region of interest and they give almost identical 
residuals. This was expected since in our region of inter¬ 
est the diffuse emission is basically accounted for by a sin¬ 
gle template based on the local Hi emission. Not surpris¬ 
ingly, the correlation functions that we obtain when using 
this model agree at the percent level with the results ob¬ 
tained with the other model, as shown in Fig. 3. For defi¬ 
niteness, we set ring_2year_P7 6_v0 . fits as our base¬ 
line model. Note that, in general, it is not recommended to 
use the model ring_2year_P7 6_vO .fits, tuned on P7 
data, with P7REP data, even though in this particular case, as 
shown in Fig. 3, the differences between the results derived 
from the two models are small. 

Finally, we have also subtracted the contributions from 
solar and lunar emission along the ecliptic. For this pur¬ 
pose we used the appropriate routines of the LAT Sci¬ 
ence Tools and selected options to obtain templates con¬ 
sistent with the data selection and IRFs choices described 
above. The model of the energy and spatial emission 
from the Sun and Moon have been taken from the re¬ 
lated papers (Abdo et al. 2011; Abdo et al. 2012), tabu¬ 
lated into the files solar_profile_v2r0 . f its^ and 

O 

lunar_profile_v2r0. fits”. 

The practical procedure we use to obtain the maps of resid¬ 
ual photon counts is to use GaRDiAn to convolve the Galac¬ 
tic emission model and the Sun and Moon templates with the 
exposure maps and the PSF which are then subtracted from 
the observed counts. Residual flux maps are then obtained by 
further dividing the residual photon counts by the exposure 
maps. 

As the Galactic diffuse emission models are not exact, 
cleaning is not perfect and the residual flux maps are still 
contaminated by spurious signal, especially on large angular 
scales. The impact on cross-correlation analyses is expected 
to be small since Galactic foreground emission is not expected 
to correlate with the extragalactic signal that we want to inves¬ 
tigate. Possible counter examples, like extinction effects, are 
small and will be discussed in Section 4 in the context of op¬ 
tical extragalactic surveys. Nonetheless we follow Xia et al. 
(2011) and apply a cleaning procedure that, using HEALPix 
tools, removes all contributions from multipoles up to £ = 10 
including the monopole and dipole. 

The resulting residual photon flux maps, which we dub £10- 
maps, for the three energy ranges considered in this work are 
shown in the three upper panels of Fig. 4. The masked area, 
also shown in the bottom plot, is represented by the uniform 
strip across the Galactic plane and further extending around 
the Fermi Bubbles and Loop I. Fluctuations have amplitude 
in the range ±1.5 • 10“^ ph cm“' s“* sr“* for the case £ > 1 
GeV according to the color code shown in the plots. Note that 

^ http://fermi.gsfc.nasa.gov/ssc/data/access/lat/ 

Model_details/FSSC_model_diffus_reprocessed_vl2.pdf 

* http://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/solar_template.html 
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Fig. 4.— Upper three plots: Fluctuations of the residual 7 -ray photon flux 
maps obtained from 60 months of Fermi-LK^ data for energies E > 500 MeV 
(top), E > \ GeV (second from the top), £ > 10 GeV (second from the bot¬ 
tom). Foreground emission from Galactic diffuse. Sun and Moon have been 
subtracted from the data as well as multipoles as large as / = 10. Different 
colors indicate fluctuations of different amplitude according to the color code 
scheme in the plots. The flat-color areas across the Galactic plane and ai'ound 
the Fermi Bubbles and Loop I correspond to the mask and have been ig¬ 
nored in the correlation analysis. The mask, which also accounts for resolved 
sources, is further shown in the bottom plot. The maps are in Galactic coor¬ 
dinates and have a resolution Vside = 512. For visualization, but not during 
the analysis, the maps have been further smoothed with 1° Gaussian filter to 
remove small scale Poisson noise. 
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Redshift z 


Fig. 5. — Redshift distributions, dN jillnz, of the different types of objects 
considered for our cross-coneiation anaiysis. SDSS DR6 QSOs (biack, con¬ 
tinuous iine), 2MASS gaiaxies (red, short-dashed), NVSS gaiaxies (magenta, 
iong-dashed), SDSS DR8 LRGs (green, short, dot-dashed) and SDSS DR8 
Main Gaiaxy Sampie (biue, iong, dot-dashed) 

for visualization, but not during the analysis, the maps have 
been smoothed with a 1° Gaussian filter to remove small-scale 
Poisson noise. The model seems to slightly over-subtract the 
7 -ray emission around (/,/?) ~ (175°,-35°), corresponding to 
the gas- and dust-rich (and thus difficult to model) Taurus 
Molecular region. Note, however, that when performing the 
cross-correlation, this region is masked by the further mask 
specific to the catalog, except in the NVSS case (see below), 
so that no bias in the results is expected from this feature. 

4. MAPS OF DISCRETE SOURCES 

In this section we describe the different catalogs of extra- 
galactic objects that we cross-correlate with the flO-maps of 
the diffuse IGRB obtained after the cleaning procedure de¬ 
scribed in Section 3. 

In this work we have considered five different catalogs: i) 
the SDSS DR 6 quasar catalog released by Richards et al. 
(2009) that should trace the ESRQ population, ii) the IR- 
selected 2 Micron All-Sky Survey (2MASS) extended source 
catalog (Jarrett et al. 2000) to trace SEGs, Hi) the NRAO VLA 
Sky Survey (NVSS) catalog of radio galaxies (Condon et al. 
1998) that we regard as alternative tracers to the ES- 
RQs, iv) the DR 8 SDSS catalog of Luminous Red Galax¬ 
ies (Abdalla et al. 2008), which should trace an intrinsically 
fainter, more local AGN population like the BE Lacs, and v) 
the DR 8 SDSS main galaxy sample (Aiharaet al. 2011) as a 
potential additional tracer of SEGs. 

The redshift distributions, dN/d\nz, of the various sources 
are shown in Fig. 5, and described in more detail in the next 
subsections. The various distributions peak at quite differ¬ 
ent redshifts, with 2MASS representing the most local pop¬ 
ulation and SDSS DR 6 QSOs the most distant one. These 
characteristics in principle enable breaking down the cross¬ 
correlation analysis in different redshift ranges, effectively al¬ 
lowing a tomographic investigation of the IGRB origin. A 
detailed description of these catalogs, except the SDSS-DR 8 
main galaxy sample, can be found in Xia et al. (201 1). Below 
we briefly summarise the main features of each sample. 


The SDSS DR 6 quasar catalog (Richards et al. (2009), 
hereafter DR 6 -QSO) contains about Vq ~ 10® quasars with 
photometric redshifts between 0.065 and 6.075, covering al¬ 
most all of the north Galactic hemisphere plus three narrow 
stripes in the south, for a total area of 8417 deg“ (correspond¬ 
ing to ~ 20% of the whole sky). The DR 6 -QSO dataset ex¬ 
tends previous similar SDSS datasets (Richards et al. 2004; 
Myers et al. 2006). The main improvements are due to the 
fact that this catalog contains QSOs at higher redshift and also 
contains putative QSOs flagged as objects with ultraviolet ex¬ 
cess (UVX objects). We refer the reader to Richards et al. 
(2009) for a detailed description of the object selection per¬ 
formed with the non-parametric Bayesian classification ker¬ 
nel density estimator (NBC-KDE) algorithm. 

In this work we used objects listed in the electronically pub¬ 
lished table with a “good object” flag in the range [0,6]. The 
higher the value, the more probable for the object to be a real 
QSO (Richards et al. 2009). We only consider the quasar can¬ 
didates selected via the UV-excess-only criteria “uvxts=r’, 
i.e., objects clearly showing a UV excess which represents 
a characteristic QSO spectral signature. After this selection 
we are left with Vq ~ 6 x 10 ® quasar candidates. 

In order to determine the actual sky coverage of the DR 6 
survey and generate the corresponding geometry mask we 
Monte Carlo sample the observed areas with a sufficiently 
large number of objects using the DR 6 database to ensure 
roughly uniform sampling on the SDSS CasJobs® website. 
Following Xia et al. (2009) we combine the pixelized mask 
geometry with a foreground mask obtained by removing all 
pixels with the g-band Galactic extinction Ag = 3.793 x E{B- 
V) > 0.18 to minimize the impact of Galactic reddening. 

The redshift distribution function dN/dz of the DR 6 -QSO 
sample in Fig. 5 is well approximated by the analytic function: 



/3 z™ 

m+1 \ _m+l 

t ( — ) Zq 



(19) 


where three free parameters values are m = 2.00, j3 = 2.20, and 
Zo = L62 (Xia et al. 2009). In addition, to calculate theoreti¬ 
cal predictions (Eq. 12) we follow Giannantonio et al. (2008); 
Xia et al. (2009) and assume a constant, linear bias model 
bs = 2.3. 


4.2. 2MASS 

The 2MASS extended source catalog is an almost-all-sky 
survey that contains ^ 770000 galaxies with mean redshift 
(z) ~ 0.072. In this work we have selected objects with 
apparent isophotal magnitude 12.0 < < 14.0, where the 

prime symbol indicates that magnitudes are corrected for 
Galactic extinction using = K 20 -A^, with A^ = 0.367 x 
E{B-V). We select objects with a uniform detection thresh¬ 
old (use-src = 1 ), remove known artefacts (cc_flag ^ a and 
cc-flag y z), and exclude regions with severe reddening, A^ > 
0.05, Schlegel et al. (1998). This procedure leaves approxi¬ 
mately 67% of the sky unmasked. The redshift distribution of 
the selected objects can be approximated with the same func¬ 
tional form used for DR 6 QSOs with parameters m = 1.90, 
fi = 1.75, and zo = 0.07 (Giannantonio et al. 2008). The value 
of the linear bias of 2MASS galaxies has been set equal to 
bs=l.A (Rassatet al. 2007). 

The possible incompleteness of the 2MASS catalog at 


4.1. SDSS DR6 QSO 


^ http://skyserver.sdss3.org/CasJobs/ 
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faint magnitudes might affect our cross-correlation analy¬ 
sis. For this reason we have repeated the analysis using 
smaller 2MASS samples with more conservative magnitude 
cuts: Xjo < 13.9, Xjg < 13.7 and Xjg < 13.5. More specif¬ 
ically, we computed the two-point cross-correlation function 
between these 2MASS maps and the Fermi E > \ GeV resid¬ 
ual map and find that the CCF results do not change signifi¬ 
cantly, i.e., the possible incompleteness of the larger catalog 
does not induce any spurious correlation signal. Therefore, in 
our analysis we use the larger 2MASS sample cut at X,g = 14. 

4.3. NVSS 

The NRAO VLA Sky Survey (NVSS, Condon et al. 1998) 
offers the most extensive sky coverage among the cata¬ 
logs considered here (82% of the sky to a completeness 
limit of about 3 mJy at 1.4 GHz) and contains 1.8 x 10^ 
sources. We include in our analysis only NVSS sources 
brighter than 10 mJy since the surface density distribution 
of fainter sources suffers from declination-dependent fluctu¬ 
ations (Blake & Wall 2002). We also exclude the Galactic 
Plane strip at \b\ < 5° where the catalog may be substantially 
affected by Galactic emissions. This additional cut is redun¬ 
dant with the one applied to the LAT maps. It is applied to 
compute the NVSS source surface density at this flux thresh¬ 
old which turns out to be 16.9 deg“^. 

The redshift distribution at this flux limit has been deter¬ 
mined by Brookes et al. (2008). Their sample, complete to 
a flux density of 7.2 mJy, contains 110 sources with 5 > 10 
mJy, of which 78 (71%) have spectroscopic redshifts, 23 have 
redshift estimates from the K-z relation for radio sources, 
and 9 were not detected in the X-band and therefore have a 
lower limit to z- We adopt the smooth parametrization of this 
redshift distribution given by de Zotti et al. (2010): 

^(z)= 1.29-1-32.37z-32.89z^-1-11.13z^-1.25z'‘. (20) 
clz 

For the bias of the NVSS sources we adopt a linear model 
bs= 1.8 (Giannantonio et al. 2012). Note that a comprehen¬ 
sive analysis of the NVSS autocorrelation function is provided 
by Xia et al. (2010). 

4.4. SDSSDR8LRG 

To sample the spatial distribution of the LRGs we use the 
photometric LRG catalog from the final imaging of SDSS 
DR8 instead of the MegaZ LRGs sample since the latter 
has an excess of power on large scales with respect to the 
ACDM model (Thomas et al. 2011). The sample used here 
was selected using the same criteria as the SDSS-III BOSS 
“CMASS” sample defined in Ross et al. (201 1). They applied 
colour cuts to account for seeing effects, dust extinction, sky 
brightness, airmass, and possible stellar contamination. 

Ho et al. (2012) further excluded regions where E(B-V) > 
0.08 for the dust extinction, when the seeing in the /-band 
> 2.0" in FWHM, and additionally masked regions affected 
by bright stars. This selection yields a sample with ~ 8 x 10^ 
LRGs and leaves approximately 22% of the sky unmasked. 

Photometric redshifts of this sample are calibrated us¬ 
ing about 100,000 BOSS spectra as a training sample for 
the photometric catalog. The resulting redshift range is 
0.45 < z < 0.65 with a mean redshift z ~ 0.5, as shown in 
Fig. 5. Also in this case, and following Ross et al. (2011); 
Hernandez-Monteagudo et al. (2014) we assume a linear bias 
parameter bs = 2.\. 



Fig. 6 . — CCFs between Fermi maps £ > 1 GeV and 2MASS galaxies 
computed using PolSpice (black diamonds) and LS (orange dots) estimators. 
Error bars are computed by using the Jackknife resampling method. 

4.5. SDSS DR8 Main Galaxy Sample 

We consider the sample of photometrically-selected “main" 
galaxies extracted from from the SDSS-DR8 catalog. The se¬ 
lection is performed according to Giannantonio et al. (2008): 
we consider objects with extinction-corrected r-band magni¬ 
tude in the range 18 < r < 21 and with redshifts in the range 
0.1 < z < 0.9. Further, we only include objects with red¬ 
shift errors smaller than = 0.5z, which leaves us with about 
4.2 X 10^ sources with redshifts distributed around a median 
value z ^ 0.35. In addition, we adopt a foreground mask to 
minimise the effect of Galactic extinction by excluding all 
galaxies within pixels in which the r-band Galactic extinc¬ 
tion > 0.18. Finally, we have about 35 million sources for 
the analyses. 

The redshift distribution dN/dz of the SDSS galaxies can 
be approximated with the same functional form used for DR6 
QSOs and 2MASS galaxies with parameters m= 1.5, /3 = 2.3, 
and zo = 0.34. Following (Giannantonio et al. 2012) we use a 
constant bias parameter bs= 1.2. 

5. CROSS-CORRELATION ANALYSIS 

In this section we describe the cross-correlation between 
the residual IGRB flux maps and the angular distribution of 
extragalactic sources in the catalogs described in Section 4. 
All the maps we use are in HEALPix format with resolution 
ofNside=512. 

Our analysis relies on the latest version v02-09-00 of Pol¬ 
Spice, a publicly available statistical tool to estimate both 
the angular two-point cross-correlation function C^®(0) and 
the cross angular power spectrum 0^(1) of any two diffuse 
datasets, f and g, pixelized on a sphere. The code is based 
on the fast Spherical Harmonic Transforms allowed by iso¬ 
latitude pixelisations and automatically corrects for the effects 
of the masks. Datasets and masks in the form of HEALPix sky 
maps are input to the code. The output consists of the angular 
two-point correlation function, the angular power spectrum 
and its covariance matrix, which account for the effect of in¬ 
complete sky coverage and from the nominal beam window 
function and average pixel function. In our calculations, we 

See http://www2.iap.fr/users/hivon/software/PolSpice/ 
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Fermi —2MASS 1 GeV 



0.1 1.0 10.0 100.0 
9 [deg] 


Fig. 7.— Comparison of the error bars of the CCF between Fermi maps 
E > \ GeV and 2MASS galaxies computed using the PolSpice covariance 
matrix and the Jackknife resampling method. CCFs for the Jackknife sub¬ 
samples are calculated with the PolSpice estimator. The thin line shows the 
unbinned CCF. 


perform the correlation analyses in the multipole and angular 
ranges f G [10,1000] and 9 G [0.1°, 100°] for CAPS and CCF, 
respectively. 

The PolSpice estimator has been described in detail and 
thoroughly tested as a tool to measure both the spectrum 
and the covariance matrix of a sky map (Szapudi et al. 2001; 
Chon et al. 2004; Efstathiou 2004; Challinor & Chon 2005). 
Since in our previous work (Xia et al. 2011) we based our 
analysis on the Landy-Szalay (LS) estimator (Tandy & Szalay 
1993) and computed errors using a Jackknife (JK) resampling 
technique (see below), we have checked the consistency of the 
two approaches and compared the outputs from the different 
estimators. 

In Fig. 6 we compare the CCF between 2MASS and Fermi 
data above 1 GeV estimated using PolSpice with the same 
CCF estimated using FS. The agreement between the two 
methods is at about 10% in the first angular bin and few % 
in the other bins, well within the amplitude of the l-cr ran¬ 
dom errors, and demonstrates that our results are robust to the 
choice of the estimator for two-point statistics. Note that, only 
in this particular case, the angular binning is linear, dictated 
by the FS estimation method which, being particularly com¬ 
putationally expensive, is applied using maps with a coarse 
Vside = 64 pixelization. 

Analogously, in Xia et al. (2011), in order to estimate the 
covariance matrix we used the Jackknife resampling method 
(Scranton et al. 2002), which divides the data into M patches 
and estimates the covariance matrix as 

M- 1 “ 
k=l 

( 21 ) 

where £,f"^{9) is the correlation function estimated in the 
k-th subsample and ^™°‘‘"(0) is the mean correlation func¬ 
tion over the M subsamples. PolSpice provides an estimate 
for the covariance matrix of the CAPS, (Efstathiou 

2004). From this, and using the procedure described in 
Planck collaboration et al. (2013) we obtain the covariance 


Fermi —2MASS CCF E>500 MeV 



x-cross 



Fig. 8 .— Comparison of Fermr-2MASS E> 500 MeV CCF (lower panel) 
and CAPS (upper panel) for two different Galactic foreground models. 

matrix for the CCF: 

,_, 7/'4-1 OP'A-i 

cll-HY. -^^P^(cos(0))P|/(cos(6»'))M^,^, , 

t V 

( 22 ) 

which is then averaged over the angular separations 9 and 9' 
within each bin to obtain a binned covariance matrix. In Fig. 7 
we show the same CCF plotted in Fig. 6 with two sets of error 
bars corresponding to the (square root of the) diagonal ele¬ 
ments of the Polspice and Jackknife covariance matrices. The 
agreement between the two sets of error bars is excellent, as 
the agreement between the off-diagonal elements (not shown) 
for which the largest difference does not exceed 10%. 

6 . VAFIDATION AND CHECKS 

In this section we assess the validity of the different steps 
of our analysis and assess the robustness of our results. For 
brevity we only present a subset of cross-correlation analyses 
involving the Fernn'-FAT and 2MASS maps. However, we 
have performed the very same robustness tests for all cross¬ 
correlation analyses described in this work. 

6.1. Test with different Galactic diffuse models 

The cleaning procedure described in Section 3 is potentially 
prone to systematic errors that may affect the correlation anal¬ 
ysis. We searched for these systematic effects by using two 
different emission models also described in Section 3 to cor¬ 
rect for Galactic emission. The results show that the corre¬ 
lation signal does not change appreciably when using either 
model. 
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Fermi —2MASS CCF E>500 MeV 




0.1 1.0 10.0 100.0 
6 [deg] 


Fig. 9.— Comparison between CCFs with and without Galactic foreground 
cleaning for the E> 500 MeV Fenni-IMASS case (upper panel) and Fermi- 
NVSS case (lower panel). 


Given the importance of this issue we performed a third test 
in which we adopted an alternative Galactic emission model 
fully based on GALPROP. In this test we used the GAL- 
PROP ‘model A’ described in the fermi-LAT collaboration 
paper (Ackermannet al. 2014a) which is one of the models 
used to assess the systematic uncertainties in the derivation 
of IGRB energy spectrum due the Galactic diffuse emission 
modeling. The model consists of two components, inverse 
Compton emission and gas emission (pion decay plus bremm- 
strahlung). Together with an isotropic template we normalize 
this 3-component model by fitting the high Galactic latitude 
7 -ray sky (more precisely using the mask of Fig. 4), leav¬ 
ing the normalization of all the 3 components free to vary 
in the fit. The tuned model is adopted as in Section 3 to 
clean the Galactic diffuse emission. The usual £10 cleaning 
procedure is then used to derive the final residual map from 
which we calculate the CCFs and CAPS. In Fig. 8 we com¬ 
pare the CCFs and CAPS between 7 rays and 2MASS for 
the gll_iem_v05_revl. fit model already used in Sec¬ 
tion 3 to the alternative model results. From the plots it can 
be seen that some difference exists at low multipoles (below 
~40). However the difference is within the 1 a error bars and 
decreases at small scales (high multipoles), where the signal 
is higher. This shows that the results of our cross-correlation 
analysis are robust with respect to the modeling of the Galac¬ 
tic diffuse emission. 

Finally, for completeness, we show in Fig. 9 the CCF be¬ 
tween the E >500 7 -ray map without any foreground clean¬ 


ing and 2MASS and NVSS. The plot shows that removing the 
Galactic emission is important to reduce the size of the error 
bars. On the other hand, even without any removal, the corre¬ 
lation is not biased, as expected given the fact that no corre¬ 
lation between the Galactic emission and LSS is in principle 
present. This in turn also implies that an imperfect Galac¬ 
tic emission removal would similarly not introduce any bias, 
although the error bars could be not optimal. 

6.2. Use of different Galactic and point sources masks 

Incorrect masking is another potential source of systematic 
errors. To tests the robustness of our results against the choice 
of the mask and to check the possible presence of Galactic 
foreground contamination we have varied the Galactic lati¬ 
tude cut from b = 20° to b = 60° in steps of Ab = 10°, as 
in Xiaetal. (2011). In addition, we performed the cross¬ 
correlation analysis using only the northern or the southern 
hemispheres of the maps. In all the cases we explored the re¬ 
sults turned out to be consistent within the errors that clearly 
increase with the size of the masked region. 

Inefficient excision of discrete sources is yet another poten¬ 
tial concern. We estimate its impact on our analysis by us¬ 
ing different masks corresponding to excluding sources from 
two different catalogs; 2FGL and a preliminary version of the 
3FGL. Moreover, to excise sources we used two criteria: i) 
we masked out circles of 1 ° radius centered on all sources 
and ii) we masked out circles of 2° radius for the 500 bright¬ 
est sources and circles of 1° radius for the others. The re¬ 
sults of the correlation analyses turned out to be insensitive to 
the choice of the source mask. Clearly, increasing the size of 
the masked areas decreases the risk of contamination but also 
decreases the significance of the correlation signal. There¬ 
fore, in an attempt to compromise between maximizing the 
statistical significance and minimizing contamination in the 
different cross-correlation analyses we proceed as follows; i) 
for the cross-correlation with the NVSS and 2MASS catalogs 
and for E > 500 MeV, which represents the case of large sur¬ 
veyed area and large contamination due to the broadening of 
the Fermi-LAJ' PSF below 1 GeV, we adopt the most con¬ 
servative source mask based on the 3FGL catalog and with 
larger (2° radius) circles, ii) for NVSS and 2MASS and £ > 1 
GeV, with a lower contamination, we use 3FGL and 1°, Hi) for 
all SDSS-based catalogs, which have the smallest sky cover¬ 
age, we have considered a less aggressive 2° 2FGL mask for 
E > 500 MeV, and iv) a 1° 2FGL mask for higher energy cuts. 

While the results are robust against the choice of the source 
mask, their significance can change appreciably. For the 
NVSS and 2MASS catalogs, using different source masks 
varies the size of the error bars by 20-30%, so the choice of 
the mask is not critical in these cases. For the various SDSS- 
based catalogs, instead, using 3FGL rather than 2FGL sig¬ 
nificantly reduces the significance of the results and, for this 
reason, we opt for the least conservative mask. 

6.3. Robustness to y-ray event conversion 

In our analysis, to maximize statistics, we consider both 
front- and back-converting 7 -ray events. However, the two 
types of events have different characteristics, most notably 
back-converting events have a larger PSF. To check whether 
the nature of the conversion affects our results we divided 
the 7 -ray datasets into two subsamples, each one contain¬ 
ing only front-converting or back-converting events. The re¬ 
sulting maps were cleaned using the same procedure, i.e.. 
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beam window functions 



Fig. 10.— Effective squared window functions of the beam, in 

the 3 energy ranges E > 500 MeV (red, dot-dashed), E > \ GeV (black, 
long-dashed), £ > 10 GeV (blue, short-dashed). 

convolving the various model templates with the appropriate 
IRFs. 

As a result of halving the number of events, the significance 
of the correlation is somewhat reduced for each subsample, in 
particular for the SDSS catalogs. Within the increased error 
bars, however, we do not observe any bias among the three 
datasets (front only, back only and front plus back). We con¬ 
clude that possible systematic differences between the front 
and back datasets are below the present statistical uncertain¬ 
ties, and we thus decided to perform the analysis using both 
type of events jointly. 

6.4. Sensitivity to the PSF of the detector 

As we shall see, a significant fraction of the cross¬ 
correlation signal observed in our analysis originates from 
small angular scales comparable with the angular extent of 
the LAT PSF. As a consequence, we need to estimate the ef¬ 
fect of the PSF and include it explicitly in our analysis. 

The PSF smears out the signal from small to large angu¬ 
lar scales, hence reducing the amplitude of both the CCF and 
CAPS at small angular separations and large multipoles, re¬ 
spectively. However this effect can be modeled if the PSF 
of the telescope, or the instrumental beam, is measured accu¬ 
rately. In the case of the LAT the beam depends on the energy, 
and the PSF can be determined either from observations by 
stacking the images of bright point sources (Ackermann et al. 
2013) or from a Monte Carlo simulation of the detector per¬ 
formance (Ackermann et al. 2012b). The characterization of 
the PSF has improved with the P7 and P7REP data release 
and a discrepancy at high energies (> few GeVs) between 
the Monte Carlo PSF and the in-flight PSF present for the P 6 
data is now significantly reduced (Ackermann et al. 2012b). 
The beam shape is part of the IRFs and can be estimated us¬ 
ing the LAT Science Tools. In particular, we used the tool 
gtpsf to obtain the PSF as a function of energy and angu¬ 
lar separation of the photon from its true arrival direction. As 
the latter is a function of one angle only we are neglecting 
the ellipticity of the beam which, in any case, turns out to be 
negligible. It is more convenient to consider the effect of the 
beam in harmonic space, where it can be expressed as a mul¬ 
tiplication rather than in configuration space, where it would 
be a convolution. Indeed, if Ci(E) represents the true CAPS at 
a given energy, then the measured one is CfE) = Wi{E)Ci{E), 


E= 0.5-100 GeV 



200 400 600 800 1000 

L 


E= 10-100 GeV 



Fig. 11.— Measured auto power spectrum of the LAT maps at |h| > 30° 
(black asterisks) and ratio between the APS and the average beam window 
squared (red, open dots) for 3 energy bands. 

where WfE) is the beam window function. The latter can be 
expressed as a Legendre transform; 

W!iE) = 2TrJ c/cos(6»)P/(cos6»)PSF(6»,£), (23) 

where Pi{x) is the Legendre polynomial of order I and 
PSF(0,£’) is the shape of the beam. Since we are analyzing 
data integrated over a fairly large energy bin within which the 
PSF can vary significantly, the effective window function for 
the bin will be a weighted average over the energy range: 

1 dN 

W,(Ei<E<E 2) =Tr dE Wi(E)—(E ), (24) 

N Je, dE 













14 


Xia et al. 


Fermi —2MASS CCF E>500 MeV 
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Fermi —2MASS CCF E>500 MeV 





Fermi-SDSS CCF E>500 MeV 


Fermi-SDSS CCF E>500 MeV 


Fermi-SDSS CCF E>500 MeV 
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Fermi-NVSS CCF E>500 MeV 


6 [deg] 

Fermi-NVSS CCF E>500 MeV 


6 [deg] 

Eermi-NVSS CCF E>500 MeV 
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Fig. 12.— Cross-correlations between the Fermi E >500 MeV map and 3 mock realization of each of the 3 catalogs 2MASS, SDSS main galaxy sample 
and NVSS, compared with the correlations with the true catalogs. The three mock realizations refer to the cases of catalog galaxies with scrambled Galactic 
coordinates {b —> —b and I —r—l) and catalog galaxies randomly distributed (MC (Monte Carlo) label in the plots) over the catalog sky-area (see text more details). 


where dN/dE represents the differential number of photon 
counts in the region of the sky we want to analyze, N = 
dE(dN/dE){E), and [Ei,E 2 ] is the energy bin considered. 
Finally, there is a further window function to take into ac¬ 
count due to the fact that we use a pixelized map, which is the 
pixel window function itself Ifpix.;. The pixel window func¬ 
tion depends on the size of the pixel (and on the shape of the 
pixel itself). It can be easily extracted using the appropriate 
HEALPix tools. The final window function is then given by 
Wi{E\<E <£ 2 ) • Wpix,/. The effective window functions for the 
3 energy ranges considered in this work are shown in Fig. 10. 

Once determined, the effective window function can be fed 
into PolSpice to recover the true CAPS and the CCF from the 
measured ones. However, we find that the algorithm used to 
perform the deconvolution is quite unstable, especially in the 
CCF reconstruction. For this reason we take the opposite ap¬ 
proach and instead of deconvolving the signal, we convolve 
the model predictions and compare them to the measurement. 
More explicitly, if C; is the model CAPS and is the es¬ 
timated effective window function in the bin AF, then the 
convolved CAPS is 


cY = CiW,^\ (25) 


and, analogously, the convolved CCF is 

CCF'^( 6 ») = Cf'PiicosB). (26) 

/ 

To test the validity of our PSF correction procedure we 
proceed as follows: i) we consider the high Galactic lati¬ 
tudes \b\> 30° region of the sky (to reduce the impact of 
the Galactic contamination) but without masking the loca¬ 
tions of known point sources; ii) we calculate the auto power 
spectrum for this region (APS, not the CAPS); Hi) we then 
apply the pipeline described in the above paragraphs to ob¬ 
tain an empirical estimate of the window function in each of 
the three energy bands considered. The resulting APS in this 
case will be dominated by the bright point sources and is ex¬ 
pected to match the Legendre transform of the window func¬ 
tion (squared): (W}^^)^. The results are shown in the three 
panels of Fig. 11. Since the APS of the map, represented 
by black asterisks, is expected to match the window function 
then the flatness of the ratio between the APS and in¬ 

dicates that our hypothesis is correct and that our estimated 
window function is robust in all three energy ranges and at all 
multipoles, apart from a small overestimation at very high 1. 
Note that, in contrast to Ackermann et al. (2012a) where the 
7 -ray APS is also considered, we neglect here the effect of 
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Ferml-2MASS CCF E>500 WeV 


Ferml-LRG CCF E>500 MeV 


Ferml-SDSS CCF E>500 MeV 


Ferml-QSO CCF E>500 MeV 


Ferml-NVSS CCF E>500 MeV 





Fig. 13.— 1-dimension profile of the l-halo-\iks normalization (in arbitrary units) from the joint l-halo-\iks and SFGsl 2-dimensional fit. Each panel 
shows the case of the fit to a single catalog and energy band CCF. 


Poisson shot noise which is sub-dominant at all multipoles, 
given the very strong APS signal from bright point sources. 
Instead, no shot noise term needs to be considered for CAPS, 
since the uncorrelated noises from the two maps being cross- 
correlated do not produce a net non-null noise CAPS, contrary 
to the APS case. The slight over-estimate of the APS PSF cor¬ 
rection in the 700-1000 1-range at the level of 20-30% turns 
into a 10-15% systematic effect for the CAPS, where the PSF 
correction is given by rather than On the other 

hand, as we will see in section 7 all CAPS are compatible 
with zero in this multipole range, except for a weak signal 
for NVSS, so that the error is dominated by statistical random 
errors. We will thus neglect the above systematic effect. 

6.5. Cosmic-ray Contamination 

The IGRB maps we obtained in Section 3 contain, be¬ 
sides true 7 -ray events, also some contamination from cos¬ 
mic rays which have been mis-classified as 7 rays. With the 
P7REP_CLEAN event selection that we used, the cosmic-ray 
contamination is at the level of 15-20% of the IGRB flux 
above 1 GeV, rising to 40-50% at 500 MeV (see Fig. 28 in 
Ackermann et al. (2012b)). Since the contaminant cosmic 
rays are not expected to correlate with cosmological struc¬ 
tures, they do not induce systematic errors in our analysis. In¬ 
stead, they will only increase random error because the signal- 
to-noise ratio of the 7 -ray signal will be reduced. Nonethe¬ 
less, to verify this hypothesis we used the IGRB maps pro¬ 
duced with the P7REP_ULTRACLEAN selection, which has 
a slightly reduced cosmic-ray contamination with respect to 
P7REP_CLEAN, at the level of 10-15% of the IGRB flux 
above 1 GeV and 30-40% at 500 MeV (Ackermann et al. 
2012b). We computed CAPS and CCFs for this case and 
found that the results are indistinguishable from those ob¬ 
tained with the P7REP_CLEAN selection. A more stringent 
test could be performed using special event selection criteria 
designed to further reduce the CR contamination for specific 
studies of the IGRB spectrum as in Abdoet al. (2010a) and 
Ackermann et al. (2014a). These selections should in princi¬ 
ple allow to reduce the error bars. However, these selection 


criteria introduce more conservative cuts to reduce the back¬ 
ground. Consequently, the benefit of the better cleaning is 
counteracted by the effective area reduction, resulting typi¬ 
cally in no effective reduction of the error bars. Indeed, the 
optimal selection should balance purity and photon statistics. 
However, the search for such compromise is beyond the scope 
of our analysis. 

6 . 6 . No-signal tests 

To check the robustness of the results we performed further 
tests using mock catalogs with no cross-correlation with LSS, 
verifying that the computed CCFs are compatible with a null 
signal. 

In Fig. 12 we show the cross-correlation between the Fenni 
E >500 MeV map and 3 mock realization of each of the 3 
catalogs 2MASS, SDSS main galaxy sample and NVSS. The 
correlations are compared with the ones with the true catalogs. 
For each catalog 2 mock realizations were built scrambling 
the Galactic coordinates of each galaxy of the sample, chang¬ 
ing b ^ -b in one case and I -I in the other. These two 
realizations preserve the intrinsic clustering of the catalog but 
remove the cross-correlation with LSS. To compute the CCF 
we use the corresponding scrambled coordinate mask of the 
given catalog. A third realization was performed creating a 
Monte Carlo catalog redistributing the galaxies of the catalog 
randomly over the sky-area covered by the catalog. In this 
case the new catalog contains no intrinsic clustering. To com¬ 
pute the CCF we use in this case the original catalog mask. 

The plots Fig. 12 show that the correlation present for the 
true catalog disappears when the mock catalogs are used, as 
expected. We note that the size of the error bars are typically 
smaller than the true catalog CCF error bars when the Monte 
Carlo catalog is used, but not in the case of the scrambled- 
coordinates catalogs. This is likely due to the fact that the 
Monte Carlo catalogs do not contain intrinsic clustering as 
opposed to the true and scrambled-coordinates catalogs, and 
emphasizes the importance of the errors cross-checks we per¬ 
formed in section 5. 

7. RESULTS 
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Fermi-QSO 500MeV 



Fermi-QSO IGeV 



Fermi-QSO lOGeV 



Fermi-QSO 500MeV 



Fermi-QSO IGeV 



Fermi-QSO lOGeV 



Fig. 14.— CAPS (upper panels) and CCF (lower panels) estimated from the SDSS DR6 QSOs map and the Fermi-hKT IGRB maps in three energy bands. The 
three panels refer to three energy cuts E > 0.5 GeV (left panels), E > \ GeV (middle panels) and £ > 10 GeV (right panels). Error bars on the data points (orange 
dots) represent the diagonal elements of the PolSpice covariance matrix. Model predictions for different types of sources are represented by continuous curves: 
FSRQs (red, dashed), BL Lacs (black, solid) star-forming galaxies (blue and green, dot-dashed) All the models ai‘e a priori models (i.e., not fitted) normalized 
assuming that the given source class contributes 100% of the IGRB. 


Fermi-2MASS SOOMeV 


Fermi-2MASS IGeV 


Fermi-2MASS lOGeV 








Fig. 15.— Analogous to Fig. 14 using 2MASS galaxies 
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In this section we show the results of our cross-correlation 
analysis both in configuration space (i.e., the angular CCF) 
and in harmonic space (i.e., the CAPS) obtained by combining 
the cleaned Ferm/-LAT IGRB maps with the angular distribu¬ 
tions of objects in the various catalogs presented in Section 4. 
The analysis is carried out in three different energy bands, and 
the results are compared with theoretical predictions obtained 
under the assumption that the extragalactic, diffuse IGRB sig¬ 
nal is generated by a combination of three different types of 
unresolved sources, BL Lacs, FSRQs and SFGs, described 
in Section 2.1. As explained in Section 6.4 we do not at¬ 
tempt to deconvolve the measured correlation function from 
the instrument PSF The predicted correlations are, instead, 
convolved with the PSF itself and then compared with the 
measurements. As we anticipated in previous sections, we 
take the relative contributions of the different types of sources 
to the IGRB as a free parameter of the model. In the next 
section we will perform a quantitative analysis to constrain 
these parameters. In this section and in the following plots, 
we use a priori models assuming that each source class con¬ 
tributes 100% of the observed IGRB spectrum, with the pur¬ 
pose of an equal-footing and quick comparison between the 
data and the different models. In fact, we assume a refer¬ 
ence “IGRB” to normalize the model predictions since, as al¬ 
ready discussed in the previous sections, the results of our 
correlation analysis do not depend on the measured IGRB 
and its uncertainties. More precisely, our reference IGRB is 
1(E) = Iq(E/Eo)-^-^ with £0 = 100 MeV and Iq = 1.44 x 10“’ 
MeV cm“^s“*sr“', which is consistent with the measured one 
(Abdo et al. 2010a). As a result, the integrated intensities of 
the IGRB in the 3 energy ranges that we consider here are 
1.0 X 10-^ 4.0 X 10“^ 1.5 X 10“* cm-^s-'sr-', for E > 500 
MeV, E > I GeV, £ > 10 GeV, respectively. 

In the following section we illustrate the results of cross- 
correlating the individual catalogs. Unlike in the analysis 
presented by Xia et al. (2011), we now observe a significant 
cross-correlation signal that, in Section 8 , we compare with 
model predictions to infer the nature of the sources that con¬ 
tribute to the IGRB. 

To asses the significance of the signals we use the usual 
likelihood ratio test assuming a gaussian likelihood L oc 
exp(-x^/ 2 ) with 

x" = - m,(/sfg))Cr ]idj - nijiU ^)), (27) 

ij 

where Cij is the covariance matrix among the different an¬ 
gular or multipole bins i computed using PolSpice, d, repre¬ 
sents the data, i.e., the CCF or CAPS measured at the bin i, 
and m,(/sfg) is the model prediction which depends from the 
parameter /sfg, i.e., the normalization of the model CCF or 
CAPS (see also Section 8 ). We use as model the SFGsl model 
with free normalization, although we note that the xlf 
the significances calculated using the other models are very 
similar. In Eq. 27 the sum extends over 10 angular bins log¬ 
arithmically spaced between 0 = 0.1° and 100° for the CCF 
and over 10 multipole bins logarithmically spaced between 
£ = 10 and i = 1000 for the CAPS. The resulting test statistics 
(TS) in this case simplifies as TS= Xo“Xbf’ where Xq is the 
of the data with respect to the null hypothesis (CCF(0)=O 
or CAPS(£)=0) and Xhf is the best-fit x^ of the data with re¬ 
spect to the model. The derived TS significances are shown in 
Table 2 for the CCFs and Table 3 for the CAPS and are com¬ 


mented in the sub-sections below for each catalog. The tables 
also report xlf and the significances in crs assuming a = '/TS. 

7.1. 1 -halo-like term 

As discussed in section 2.3, a further contribution to the 
cross-correlation can arise from a 1-halo term or if part of 
the sources of a given catalog are also themselves 7 -ray emit¬ 
ters. We denoted these terms collectively as 1 -halo-\\ke. To 
test empirically the possibility of the presence of a f-ZraZo-like 
term we adopt the following procedure. For each catalog and 
for each energy band we perform a two-parameter fit using 
a similar x^ as in Eq. 27, but modeling m,- as the sum of the 
SEGsl model with free normalization plus a further term pro¬ 
portional to the PSE profile, i.e., oc for the CAPS, and to 
the related harmonic transform for the CCF The latter is rep¬ 
resentative of a correlation at zero angular separation which is 
spread at larger angular scales by the effect of the PSF, as ex¬ 
pected for a f-ZzaZo-like term. Again, the results change only 
marginally if a model different from the SFGsl model is used. 

In Fig. 13 we show the 1-dimensional x^ of the normal¬ 
ization of the 7-ZraZo-like term profiled over the remaining 
parameter, i.e., the function obtained from X^(/sfgi/ih) after 
minimizing over /sfg. The plots refer to the fit to the CCFs. A 
fit to the CAPS gives similar results. It can be seen that for 
all the energy bands and the four catalogs SDSS DR 6 QSO, 
2MASS, SDSS LRGs and SDSS MG the significance of this 
extra term is typically below 1 a, reaching the largest signifi¬ 
cance of little more than 1 a only in the case of the correlation 
of 2MASS with 7 rays above 1 GeV. The only exception is the 
NVSS case, where, clearly, a strong preference for a 1-halo- 
like term is present. We will discuss further the NVSS case in 
the dedicated section below, while, given the lack of signifi¬ 
cant indications of the presence of 7-ZiaZo-like contributions, 
we will not consider the other catalogs further in the following 
and in the global fitting described in Section 8 . 

7.2. Cross-correlation with SDSS DR6 QSOs 

The DR 6 QSO catalog contains AGN at high redshifts that 
should preferentially trace bright FSRQs, whose redshift dis¬ 
tribution also extends to high redshifts. The results of the 
cross-correlation analysis are shown in Fig. 14. For readabil¬ 
ity, in all the plots of this and the following sections, we only 
show the predictions for the BLLacsl model, since the pre¬ 
dictions for the BLLacs2 model are rather similar. The up¬ 
per panels show the CAPS in three energy bands (increasing 
from left to right). We plot the corresponding CCFs in the 
lower panels. At small angular scales 0 < 1° we observe a 
cross-correlation signal which is more significant in the low 
energy band (^ 4.5 a for the CCF and ~ 5 cr for the CAPS), 
where photon statistics are higher. The fact that a weaker sig¬ 
nal (2-3 a) is also present for £ > 1 GeV suggests that the 
cross-correlation is genuine and not an artifact from system¬ 
atic errors in the cleaning procedure. 

The observed CCF is perfectly consistent with the theoret¬ 
ical predictions of all the a priori models considered: BL 
Lacs, FSRQs and SFGs for all three energy ranges. This is 
not entirely surprising since the dl(> E)/dz of all these mod¬ 
els overlap significantly with the dN/dz of the DR 6 QSOs. 
The similarity of the model predictions implies that BL Lacs, 
FSRQs, SFGs and DR 6 QSOs all trace the underlying mass 
density field at high redshifts. 

We note that the SDSS DR 6 catalog of QSOs is prone to 
a systematic error that we did not investigate in the previous 
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Fig. 16.— Analogous to fig. 14 using NVSS galaxies 
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Fig. 17.— Analogous to fig. 14 using SDSS DR8 Luminous Red galaxies 


sections: contamination by stars. To investigate this issue and 
assess the magnitude of the effect, we have extracted a large 
number 8 x 10'*) of SDSS DR6 stars with apparent mag¬ 
nitudes in the range 16.9 <g< 17.1 from the CasJobs web¬ 
site. We then estimated the cross-correlation between this star 
catalog and the Fermi maps. The resulting CCFs turned out 
to be consistent with zero, showing that any residual stellar 
contamination does not correlate with the IGRB and does not 
contribute to the cross-correlation signal. 


7.3. Cross-correlation with 2MASS galaxies 

The 2MASS survey catalog is the most local sample that 
we have considered. These near-infrared-selected galaxies are 
likely to trace the local SFG population rather than the AGN 
population. The results of the analysis are shown in Fig. 15. 
We observe a signal in the CCF at angular separations smaller 
than ^10° with a significance ^ 3.5(7, and a signal in the 
CAPS with a similar 3.5 a significance which appears to re¬ 
sult mainly from multipoles smaller than ~ 200. The angu- 
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lar extent and the amplitude of this signature depends on the 
energy band. Intriguingly the significance remain stable at 
£ > 1 GeV, even slightly increasing, especially in the case of 
the CCF, possibly indicating a signal peaking at around GeV, 
as expected for the case of nearby SFGs. 

The comparison with the models excludes, with high signif¬ 
icance, that BL Lac could give a dominant contribution to the 
IGRB diffuse emission at low redshift. In this respect both the 
CCF and the CAPS provide strong constraints. This result is 
in agreement with similar independent findings, based on the 
population studies of resolved BL Lacs (Abdoet al. 2010b) 
and the anisotropy of the IGRB (Cuocoetal. 2012), which 
both indicate a low contribution BL Lacs (<20-30%) at least 
up to ~50 GeV. Above 50 GeV the contribution is more un¬ 
certain and can be more significant (Ajello et al. 2014). In 
this respect, our constraints in the energy range £ > 10 GeV 
are weak and do not provide a direct test. Conversely, both 
SFGsl and FSRQs are equally good candidates for the IGRB 
in the energy range explored here. The close match with the 
data stems from the fact that they have similar dl{> E)/dz 
at low redshifts. As expected, the predictions for the SFGs2 
model are significantly different from the SFGsl one and do 
not fit the data. This implies that in model SFGs2 the con¬ 
tribution from star-forming galaxies should be <20-30% of 
the total, as for the BL Lacs. As discussed in Section 2.1 
the large differences between the two SFG models originates 
from the different distributions dl(> E)/dz with the distribu¬ 
tion for SFGs2 peaking at very low redshift as opposed to the 
SFGsl one that extends to high redshift. The implications of 
these differences are discussed in Sections 8 and 9. 

7.4. Cross-correlation with NVSS galaxies 

Fig. 16 shows the results obtained by cross-correlating 
Fermi maps with the NVSS galaxy catalog. In this case we de¬ 
tect a CCF signal with a strong significance of about ~ S.Ocr 
both for £ > 500 Me V and £ > 1 GeV on small (6 < 1 °) angu¬ 
lar scales. In fact, we detect a strong signal also in the highest 
energy bin (^5.0cr), though only at 0 < 0.2°. The fact that the 
peak in the CCF narrows with increasing energy is quite in¬ 
formative and indicates that the signal is intrinsically confined 
to very small 9 and it extends to larger 6 values only because 
of the spreading out effect by the LAT PSF, especially at low 
energies. The width of the peak, ^ 1.5°, ~ 1.0° and ~ 0.2° 
for the CCF at £ > 0.5,1,10 GeV is, indeed, also compatible 
with the width of the LAT PSF at these energies. The CAPS 
gives similar significances and provides additional informa¬ 
tion on the characteristics and possible origin of the signal. 
Differently from the CAPS with other catalogs, in fact, the 
CAPS with NVSS is characterized by a strong signal at very 
high multipoles (up to / ^ 1000 ) confirming that the signal 
comes mostly from small scales. 

All models provide a good match to the data at large {6 > 
1°) angular scales. At smaller separations, however, the ob¬ 
served signal overshoots model predictions, especially in the 
SFGsl case, as confirmed by the high Xm in Tables 2 and 
3. This excess signal correlation on small scales does not 
seem to be related to the large-scale clustering of astrophys- 
ical sources. Instead, this correlation seems to be better de¬ 
scribed by a J-halo-like term. Indeed, as seen in section 7.1, 
the NVSS case is the only one where a J-halo-like term is 
strongly detected. Using a J-halo-like term as an alternative 
model to calculate the significance of the signal improves the 
quality of the fit to both CAPS and CCF, as confirmed by the 


decrease in the Xbf values and the corresponding increase in 
statistical significance to ^ 10.0 cr. 

It is unclear if this small-scale signal is due to a pure 1- 
halo term or to the possibility that a significant fraction of 
NVSS sources might also be 7 -ray emitters. Indeed, the fact 
that NVSS sources are known to be good candidates for 7 - 
ray emission and that this catalog is routinely searched to 
identify the counterparts of 7 -ray sources (Nolan et al. 2012; 
Acero et al. 2015) is an argument in favor of the second pos¬ 
sibility. 

A further explanation, possibly not entirely independent 
of the previous ones, is the presence of duplicate objects 
in the NVSS catalog. It is well known that a large frac¬ 
tion of close pairs are in fact single objects with a promi¬ 
nent radio jet wrongly classified as a separate, companion 
object (Overzier et al. 2003). The net result is an excess of 
pairs at small angular separations which is responsible for 
an unphysical large auto-correlation signal at small angles 
(Overzier et al. 2003), and could thus induce a corresponding 
cross-correlation excess. 

For all the above reasons, we adopt a conservative approach 
and consider the NVSS cross-correlation at angles 0 < 1° and 
multipoles I > 100 as arising from physical processes that are 
not associated to the large-scale structures and will ignore it 
in the x^ analysis performed in the next section. 

7.5. Cross-correlation with SDSS DR8 LRGs galaxies 

The results of the cross-correlation analysis between the 
SDSS DR 8 LRGs and the Fermi-LPJT maps are shown in 
Fig. 17. In this case we do not detect any correlation sig¬ 
nal. In fact the CCF drops below zero at very small angu¬ 
lar separations, although the significance of this feature is 
very weak. A possible reason for this surprising behaviour 
is the aggressive procedure used to remove possible system- 
atics from the raw LRGs data (Hoet al. 2012) which might 
remove the genuine correlation signal together with the spuri¬ 
ous one. At larger (9 > 0.2°) separations the correlation signal 
is consistent with zero. This is in agreement with the model 
predictions for SGFsl and, to a lesser extent, to BLLacsl and 
SFGs2. This is not surprising since the dl{> E)/dz of these 
sources barely overlap with the narrow dN jdz distribution of 
the LRGs. On the contrary, the FSRQ model predicts a sig¬ 
nificant cross-correlation, which is at variance with the data. 

7.6. Cross-correlation with SDSS-DR8 main galaxy sample 

In Fig. 18 we plot the estimated CAPS and CCF between 
the £emu-LAT maps and the SDSS-DR 8 galaxies in the main 
sample. We observe a correlation signal at small angles at 
about > 3 (T level for both the £ > 500 MeV and the £ > 1 
GeV cases, similarly to the 2MASS case. The observed CCF 
is marginally consistent, for £ > 1 GeV, with theoretical pre¬ 
dictions if the sources of the IGRB are SFGs in model SFGsl. 
In all AGN-based models the predicted cross-correlation sig¬ 
nal is much higher than the observed one. This is similar to 
the 2MASS case except that the dN/dz of the DR 8 galaxies 
peaks at significantly higher redshift than 2MASS galaxies. 
We conclude that in this case SFGs provide a significant con¬ 
tribution to the IGRB not only locally but also at z ~ 0.3 and 
that their contribution is more important than that of BL Lacs 
and FSRQs. In the case of the SFGs2 model, instead, SFGs 
are predicted to have a small contribution, similar to the one 
of blazars. 

8 . ANALYSIS 
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Fig. 18.— Analogous to fig. 14 using SDSS DR8 main galaxy sample 


TABLE 2 

Significance of the CCFs cross-correlations for each energy bin and catalog calculated using the SFGs 1 model with free 

NORMALIZATION. FOR EACH CASE, THE BEST FIT THE SIGNIFICANCE cr AND THE TEST STATISTICS TS VALUES ARE REPORTED. EACH FIT HAS 9 
DEGREES OF FREEDOM (10 BINS - 1 FREE PARAMETER). FOR THE NVSS CASE A FURTHER MODEL, PSF, IS TESTED. 


CCF 

2MASS 

SDSS-MG 

SDSS-LRG 

SDSS-QSO 

NVSS (LSS) 

NVSS (PSF) 



a 

TS 


cr 

TS 


a 

TS 

Xm 

cr 

TS 

xii 

cr 

TS 


cr 

TS 

E > 500 MeV 

6.2 

3.6 

12.9 

2.6 

2.7 

7.4 

4.5 

0.3 

0.1 

9.0 

4.5 

21 

30.2 

8.0 

64.9 

3.6 

9.9 

97.3 

E > 1 GeV 

10.6 

4.4 

19.4 

2.1 

3.0 

9.3 

4.6 

0.4 

0.2 

3.5 

2.3 

5.1 

45.1 

8.6 

73.6 

4.9 

10.3 

106.4 

£ > 10 GeV 

2.0 

2.1 

4.5 

6.2 

0.7 

0.5 

2.6 

0.2 

0.1 

4.8 

1.6 

2.6 

40.4 

5.1 

25.6 

5.8 

7.7 

59.4 


TABLE 3 

Same as Table 2 but using CAPS. 


CAPS 

2MASS 

SDSS-MG 

SDSS-LRG 

SDSS-QSO 

NVSS (LSS) 

NVSS (PSF) 


Xbf 

cr 

TS 

Xhf 

cr 

TS 

Xhf 

cr 

TS 

Xhf 

cr 

TS 

Xhf 

cr 

TS 

Xilf 

cr 

TS 

E > 500 MeV 

8.3 

3.4 

11.5 

4.5 

3.5 

12.1 

3.5 

0.0 

0.0 

9.7 

5.3 

28.6 

30.1 

8.3 

71.3 

7.3 

9.6 

92.3 

E > 1 GeV 

3.7 

3.6 

12.8 

3.9 

3.3 

11.2 

5.4 

0.4 

0.2 

7.6 

3.3 

10.9 

23.1 

8.4 

70.7 

5.3 

9.1 

82.8 

£ > 10 GeV 

5.1 

1.6 

2.7 

8.4 

0.7 

0.6 

4.4 

0.7 

0.5 

4.6 

2.7 

7.3 

21.0 

3.4 

11.8 

9.3 

4.8 

23.2 


To quantify the qualitative conclusions drawn from the in¬ 
spection of the correlation analysis performed in the previous 
section we now perform a comparison between model pre¬ 
dictions discussed in Section 2 and the CCF and CAPS esti¬ 
mates presented in Section 7. The aim is to estimate the free 
parameters of the models, i.e., to quantify the relative con¬ 
tribution of different types of potential sources to the IGRB 
and to assess the goodness of the fit, from which we can infer 
which is the most likely mix of source candidates responsi¬ 
ble for the observed IGRB. Here we present only the results 
of the CCF analysis since those obtained with the CAPS are 
fully consistent with those shown below. 

For each CCF estimated by comparing a galaxy catalog and 
a Fermi-LAT map above a given energy threshold we compute 


the following statistics; 

~ J ~ ™’ ( 28 ) 

ij 

where Co.g. is the covariance matrix computed using PolSpice 
that quantifies the covariance among different angular bins 6*,, 
di represents the data, i.e., the CCF measured at the angu¬ 
lar bin i, and OT,(a) is the model prediction which depends 
from a set of parameters a. We note that it is important to 
use the full covariance matrix since the different bins are sig¬ 
nificantly correlated, a feature which is typical of CCF mea¬ 
surements. Instead, the covariance matrix of the CAPS is to a 
better approximation diagonal (although some sizable corre¬ 
lations are nonetheless present, in particular for low and high 
multipoles), at the price, however, of making the interpreta- 
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TABLE 4 

Minimum for the one-, two-, and three source models and best fit values for the free parameters corresponding to the 

FRACTION OF THE IGRB CONTRIBUTED BY SFGS, BL LACS AND FSRQS. 



BLLacsl 

BLLacsl 



/sFGs 

/blls 

/fSROs 

T 

X" 

/sFGs 

/bLLs 

/fSROs 

SFGsl 

35.3 

0.72 



35.3 

0.72 



BLLacs 

44.3 


0.08 


43.1 


0.18 


FSRQs 

48.8 



0.24 

48.8 



0.24 

SFGls + BLLacs 

35.3 

0.72 

0.0 


35.3 

0.72 

0.0 


FSRQs + SFGsl 

35.3 

0.72 


0.0 

35.3 

0.72 


0.0 

FSRQs + BLLacs 

42.0 


0.06 

0.10 

43.1 


0.18 

0.0 

FSRQs + BLLacs + SFGsl 

35.3 

0.72 

0.0 

0.0 

35.3 

0.72 

0.0 

0.0 

FSRQs + BLLacs + SFGs2 

41.7 

0.14 

0.0 

0.12 

41.7 

0.14 

0.0 

0.06 


tion less intuitive since it lacks the immediate identification 
of the scale(s) responsible for the correlation, which is instead 
present for the analysis in real space with the CCFs. Thanks 
to the fact that we are considering 7 -ray flux maps rather than 
fluctuation maps we can express the model CCF as a sum of 
different contributions corresponding to the CCF of different 
source types: m,- = /^fg Csfg(6»;) -F fuuc CuudSi ) + /tsrq Cfsrq(6»;), 
where cdOi ) is the model CCF for a given type of source 
when it represents 100 % of the IGRB and fa is a free pa¬ 
rameter that quantifies the actual IGRB fraction contributed 
by the source. Note that, in our analysis, we do not require 
that X] /a = 1 ■ Instead we verify that this condition is verified 
a posteriori. 

In Eq. 28 the sum extends over 10 angular bins logarithmi¬ 
cally spaced between 9 = 0.1° and 100°. We use logarithmic 
bins to emphasize small scales where the signal-to-noise is 
higher while the choice of 10 bins is dictated by the compro¬ 
mise between the need to robustly invert the covariance ma¬ 
trix (an operation that becomes unstable when too many cor¬ 
related bins are considered) and that of maximizing the avail¬ 
able information. The number of bins used in the yf analysis 
(10) is smaller than that used in the CCF plots shown in the 
previous Section (20), in which the aim was to illustrate the 
qualitative agreement between models and data. 

The total accounts for the contributions from the CCF of 
all catalogs. The only exception is the CCF of NVSS sources 
for which we have ignored separations 0 <\° since, as dis¬ 
cussed in the previous section, the signal in that range is likely 
not related to large-scale structure clustering. In particular, the 
total for the E >500 MeV band is the quantity that we use 
to perform the bulk of the analysis detailed below. However, 
we have also considered the cases in which the x^ only ac¬ 
counts for the CCF of a subset of catalogs and/or of different 
energy bands. Investigating the contribution to x^ by differ¬ 
ent catalogs is important to illustrate the tomographic nature 
of our analysis. On the contrary, considering different energy 
cuts turns out to be not very revealing. In principle, breaking 
out the total x^ by energy bands provides additional infor¬ 
mation to identify the contribution to the IGRB by different 
sources. However, the constraints derived from the 2 higher 
energy bands are weak and the x^ discriminating power is 
dominated by the E > 500 MeV band. The practical outcome 
is that the results obtained by considering photons with £ > 1 
GeV or £ > 10 GeV are fully consistent with those obtained 
with the £ > 500 MeV energy band. 

We have performed our x^ analysis in three steps in which 
we increase the complexity of the IGRB model: i) the one- 


source scenario, in which we assume that only one type of 
source, FSRQs, BL Lacs, or SLGs, contributes to the IGRB, 
ii) the two-source scenario, in which we allow for two possi¬ 
ble contributors to the IGRB and Hi) the three-source scenario 
in which LSRQs, BL Lacs, and SLGs can contribute to the dif¬ 
fuse background. As we already noted in all three cases the 
overall normalization is free, i.e., we do not impose that the 
overall contribution should sum up to, or not exceed the ob¬ 
served IGRB. Instead we have checked that, after minimizing 
the x^, this condition is satisfied in all cases explored. 

The results of the x^ minimization are summarized in Table 
4 for the one-source (upper part), two-source (middle section) 
and three-source (bottom) scenarios. In the Table we list the 
minimum x^ value together with the three best fit parameters, 
in parentheses, i.e., the IGRB fraction contributed by SLGs, 
BL Lacs and LSRQs, respectively. In the one- and two-source 
scenarios the values of the sources not considered in the fit are 
set equal to zero, and the related space in the table is left blank 
for clarity. The two columns refer to the two different BL 
Lac models that we have considered (BLLacsl and BLLacs2). 
Note that we quote total x^ values rather reduced ones, since 
it is not straightforward to calculate the number of degrees of 
freedom involved. This quantity, in fact, is not simply equal to 
the number of bins over which the x^ is calculated due to the 
presence of correlation among different catalogs, since their 
redshift distributions and angular coverages overlap signifi¬ 
cantly. Instead, to assess the goodness of fit, we quote, for 
the case of the three-parameter models, the best-fit x^ values 
for the cross-correlation between each single catalog and the 
Fermi E >500 MeV map. These x^ values are presented in 
Table 5, which can be compared with the number of degrees 
of freedom, given approximately by the number of bins used 
minus the number of fit parameters. The results indicate that 
the fit to each catalog for the three-parameters models is ad¬ 
equate except for a tension with the QSO CCL in the SLGsl 
model that is even more prominent in the SLGs2 model. The 
tension among the models results in the under-prediction of 
the observed correlation. 

The main results of the x^ analysis are: 

• All models including a contribution from SLGs provide 
a significantly better fit than those in which the IGRB 
is contributed to by AGN only. 

• Model SLGsl performs better than SLGs2. The main 
issue with the SLGs2 model is that it provides a poor 
fit to the CCL of the SDSS-QSOs, as indicated in Table 
5, while for all other datasets the two SLG models fit 
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Fig. 19.— Plot matrix showing the 1-dimensional profile likelihoods for each component and contours of the 2-dimensional profile likelihoods for the three- 
component fit (BL Lacs, FSRQs and SFGs) to all the experimental CCFs (i.e., all catalogs and all energy ranges). The plots refer to the model BLLacsl and 
SFGsl. The plot in the upper-right corner shows the profile likelihood for the total IGRB fraction. 


TABLES 

Contribution to the best fit from the single catalog CCFs with the E >500 MeV 7-ray map, for the two models FSRQs + 
BLLacs 1 + SFGs 1 AND FSRQs + BLLacs 1 + SFGs2 . The numbers in parenthesis are the number of 9 bins used to calculate the x^ ■ 



2MASS (10) 

NVSS (6) 

SDSS-MG (10) 

SDSS-LRG (10) 

SDSS-QSO(IO) 

FSRQs + BLLacsl + SFGsl 

6.4 

1.5 

3.6 

7.7 

16.1 

FSRQs + BLLacsl + SFGs2 

6.2 

1.5 

3.1 

6.6 

24.3 


the data equally well, although with different overall 
normalizations of the SFG signal. 

• In all models explored the IGRB contribution from 
AGN is subdominant. When SFGs are included among 
the IGRB sources, the AGN contribution is consistent 
with zero. The consistency with zero simply reflects 
the limited accuracy of our analysis which does not ac¬ 
count for the fact that, based on the observed number 
count distribution of the resolved 7 -ray sources, some 
contribution from AGN is to be expected (Ajello et al. 
2012, 2014). 

• BLLacsl and BLLacs2 models have similar values 
although the normalization of the BLLacs2 component 
is approximately a factor of 2 higher than the BLLacsl 
model. 


The uncertainties in the estimates of the parameters can be 
appreciated from the sets of panels shown in Figs. 19 (SFGsl 
model) and 20 (SFGs2 model). We do not show results for 
the BLLacs2 case, since they are very similar to those of the 
BLLacsl case, when one rescales the BL Lac component by a 
factor of ~2. Among the plots, those with the 1-dimensional 
represent the contribution to the IGRB from a specific type 
of source that we obtain profiling (Rolke et al. 2005) over the 
other contributors. For example in the case of SFGs, this is 
the function obtained after minimizing the x^(/sfgj/biiac,/fsrq) 
with respect to /biiac and /fsrq. The plots show the x^ together 
with the 2 tr significance level. The plot in the upper-right 
corner also shows the derived quantity /tot, i.e., the total IGRB 
fraction, /^t = /sfg +/biiac +/fsrq- The 2D contours refer to the 
function obtained by profiling over only one parameter. In 
this case the contours are drawn in correspondence to the 1 , 
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Fig. 20.— Same as Fig. 19 but for models BLLacsl and SFGs2. 

2 and 3 a confidence levels. The constraints on the SFGs are 
rather broad; they show that, within 2 tr, the contribution to 
the IGRB of these sources varies between 25-95 %. The con¬ 
straints are tighter in the case of BLLacs2 model (50-95 %) 
not shown in the plots. A scenario with no SFG contribution 
is rejected with high statistical significance. The contribution 
of AGNs is consistent with zero but, within 2 cr, can be as 
large as 5-7% and 7-8% for the FSRQs and BL Lacs, respec¬ 
tively. 

For the SGFs2 model, which provides a worse fit than 
SFGs 1, we only obtain upper limits; the SFG contributes < 20 
% to the IGRB and the contributions of BLLacsl and FSRQs 
are limited to 10% and 20%, respectively. Note that in this 
case the contribution from FSRQs is larger than in the SFGsl 
model. Another difference between model SFGsl and SFGs2 
is the total contribution to the IGRB. The three-source model 
SFGsIh- BLLacsl + FSRQs model is able to account for a 
large fraction of the IGRB, about 70-80%, while the model 
SFGs2h- BLLacsl + FSRQs model can only be responsible 
for ~20-30% of the IGRB. Table 5 emphasizes the main is¬ 
sue with the model SFG2 h- BLLacsl + FSRQs, which is the 
poor fit to the QSO CCF Indeed, visual inspection of the CCF 
confirms that this model cannot explain the amplitude of the 
measured cross-correlation signal. Instead, model SFGsl-t 
BLLacs 1 + FSRQs provides a better fit to the data apart still a 
small residual underestimate of the QSO correlation signal. 

Finally, another interesting feature of the 2D contours 


is the non-negligible correlation among the contributing frac¬ 
tions. The fact that these contours are not completely de¬ 
generate, however, is a non-trivial result that we have ob¬ 
tained by cross-correlating the Fermi maps with different 
catalogs of objects spanning different redshift ranges or, in 
other words, to the tomographic nature of our analy¬ 
sis. To illustrate this point, let us consider the simple two- 
source model BLLacs2H-SFGsl instead of the 3-source one, 
BLLacs2H-SFGslH-FSRQs. The advantage is that in this case 
the 2D function x^(/sfgj/biiac) encodes all information which, 
instead, is partially lost when one profiles the three-source 
model. We show in Fig. 21 the 1 and 2 a contours of the 
SFGh-BL Lac contributions superimposed to the 1 and 2 cr 
contours obtained when only one type of catalog is consid¬ 
ered. Here we show the y^ values obtained by cross corre¬ 
lating the Fermi maps with SDSS QSOs, SDSS galaxies and 
2MASS galaxies. Individual constraints are fully degenerate, 
as one can only constrain the ratio of the two contributions. 
In particular, constraints from objects at low redshifts, like 
2MASS and SDSS galaxies, would only narrow the width of 
the uncertainty strip, without removing the degeneracy. It is 
only when we consider SDSS QSOs which convey informa¬ 
tion on clustering at high redshifts that we are able to remove 
part of the degeneracy. Note that the combined constraints 
are consistent with those obtained from the analyses of the 
individual catalogs. 
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shift (z > 1) clustering that is traced by the QSOs in the 
SDSS catalog. 



BLL fraction 

Fig. 21.— contours for the 2 component model SFG+BLLacs2 with re¬ 
spect to the 3 subsets of CCFs with the SDSS QSOs (pale brown, dot-dashed), 
the 2MASS (red, solid) and the SDSS MG (pink, dashed) catalogs, and the 
combined dataset (blue, solid). The complementary tomographic informa¬ 
tion from the different catalogs helps to break the degeneracies present when 
using a single catalog. 

9. DISCUSSION AND CONCLUSIONS 

In this work we have cross-correlated the Fermi-hP^ sky 
maps obtained in 60 months of observations with the angular 
positions of several types of extragalactic objects at different 
redshifts. The aim is to constrain the origin of the IGRB un¬ 
der the hypothesis that it is constituted by unresolved astro- 
physical sources that can be traced by but do not necessarily 
coincide with the objects in the catalogs. The benefit of per¬ 
forming a cross-correlation analysis of the IGRB rather than 
considering its mean amplitude or its auto-correlation proper¬ 
ties is that the cross-correlation is less prone to systematic ef¬ 
fects that may arise from the inaccurate cleaning of the 7 -ray 
maps, such as imperfect subtractions of the diffuse Galactic 
foreground, contributions from charged particles and so on. 

For this purpose we rely on two complementary statistical 
tools; the angular two-point correlation function and its Leg¬ 
endre transform: the angular power spectrum. The results of 
our cross-correlation analysis were compared with theoretical 
predictions in which one assumes that the IGRB is consti¬ 
tuted, in full or in part, by any of these potential candidates; 
SFGs, BL Lacs and FSRQs. 

The main results of our analysis are: 

• We observe a significant (> 3-4 cr) signal in the angular 
cross-correlation function of 2MASS galaxies, NVSS 
galaxies, and QSO with the IGRB on scales smaller 
than 1°. A weaker signal, ^ 3cr, is also observed for 
SDSS main sample galaxies. While in the case of 
2MASS the cross-correlation signal is observed in all 
energy bands and seems to be genuinely related to the 
underlying clustering properties of matter, in the case of 
NVSS we interpret the CCF signature as not originating 
from the large-scale structures. The NVSS signature is 
likely attributed, at least in part, to undetected 7 -ray 
sources that have counterparts in the NVSS catalog and 
to spurious close pairs in the catalog that are, in fact, 
a single object. The fact that a cross-correlation sig¬ 
nal on small scales is also observed when we consider 
SDSS galaxies and SDSS QSOs is a very interesting 
result as it suggests that the CCF signal does not solely 
originate at redshift z < 0.1 , where 2MASS and SDSS 
galaxies are found, but is also contributed by high red- 


• The fact that we now observe a signal in several cross¬ 
correlation analyses is beyond the original expectations 
of Xia et al. (2011) who performed an analysis similar 
to the one presented here using Fernu'-LAT two-year 
maps. Their forecast, however, was based on the as¬ 
sumptions that errors were dominated by Poisson noise 
from discrete photon counts. Our positive result relies 
on several improvements that have enabled us to effi¬ 
ciently remove potential sources of systematic and ran¬ 
dom errors. In particular: i) in the map cleaning proce¬ 
dure we have used three different models for the Galac¬ 
tic diffuse foreground that update and improve the one 
used in the original analysis, ii) we were able to excise a 
larger number of individually resolved sources from the 
7 -ray maps thanks to the most recent LAT source cata¬ 
logs, Hi) the addition of another set of discrete sources, 
the SDSS Main Galaxy catalog, in our cross-correlation 
analyses and iv) a better characterization of the PSF of 
the LAT. The latter improvement is probably the most 
crucial since it not only allows us a better comparison 
between model and data, but also allows to push our 
analysis down to 500 MeV, significantly increasing the 
photon statistics and reducing the amplitude of statisti¬ 
cal errors, and to smaller angular scales than Xia et al. 
( 2011 ), below 1 °, where the signal is the most promi¬ 
nent. 

• We have verified with a series of dedicated tests that 
the results of our cross-correlation analysis are robust 
to i) the cleaning procedure of the Fermi maps, ii) the 
subtraction of the Galactic Diffuse Foreground, Hi) the 
removal of the resolved 7 -ray sources, iv) the choice 
of the mask, v) the 7 -ray conversion layer in the LAT, 
vi) the statistical estimator used to measure the angular 
cross-correlation function and vH) the method adopted 
to assess the uncertainties in the CCF and CAPS and 
their covariance. In addition, we have verified that our 
characterisation and treatment of the PSF of the tele¬ 
scope is good and does not introduce any significant 
systematic error in the comparison between models and 
data. 

• The comparison between measured cross-correlation 
signals and model predictions indicates that the best- 
fit to the data is obtained when SFGs are the main, if 
not the only, contributors to the IGRB (possibly de¬ 
generate with MAGN, see below) and AGN provide a 
minor, possibly negligible, contribution. We have ex¬ 
plored different combinations of sources and different 
models for the 7 -ray contribution from BL Lacs and 
SFGs. Models that include SFGs outperform those that 
consider AGN only. And among the SFG models ex¬ 
plored, the one proposed by Fields et al. (2010) that in¬ 
cludes the effect of gas quenching and its redshift de¬ 
pendence provides a better fit to the data than the one 
proposed by Tamborra et al. (2014) which, instead, ig¬ 
nores this effect. 

Our analysis makes these statements more quanti¬ 
tative and shows that, for the model that provides the 
best-fit, SFGs contribute to 72^3^ % (but see our dis¬ 
cussion of MAGN below) of the IGRB (2 cr confidence 
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interval) whereas BL Lac and FSRQs provide similar 
contributions ranging from 0 to 8 % each. In none of our 
best fit models does the contribution to the IGRB total 
up to 100 %. This is an interesting result that keeps 
open the possibility that other types of sources could 
contribute to the 7 -ray background. In the framework 
of the cross-correlation analysis one can only speculate 
on the nature of these sources. Among the different 
options, the possibility that they consist of astrophysi- 
cal sources at high redshifts, that would not be detected 
by our cross-correlation analysis, or that they originate 
from the annihilation or decay of dark matter particles, 
are especially intriguing and will be investigated in fu¬ 
ture analyses. 

Model predictions depend on a number of parameters, 
including the bias relation of the mass tracer. Cur¬ 
rent models of galaxy evolution do not provide reli¬ 
able predictions for the bias of BL Lacs, FSRQs and 
SFGs which are only weakly constrained by observa¬ 
tions. For this reason we run a series of robustness tests 
in which we have considered the alternative bias mod¬ 
els described in Section 2.1. 

For BL Lacs and FSRQs we have considered the case 
of constant bias bfSRQ = bsLLac = 1-04 as well as that 
of a z-dependent bias matching that of a IO'^Mq dark 
matter halo. The first scenario predicts a larger cross¬ 
correlation signal for low-redshift objects (i.e., 2MASS 
and SDSS galaxies) than in the reference case. It is a 
~ 20% effect that improves the match between FSRQ 
model and data. At higher redshift the cross-correlation 
slightly decreases. However the effect is very small 
and doesn’t affect the outcome of the model vs. data 
comparison. In the second scenario the bias is system¬ 
atically larger than the reference one at all redshifts, 
significantly increasing the amplitude of the predicted 
cross-correlation. The net result is that, in this rather 
extreme case, our conclusion that SFGs contribute to 
the bulk of the IGRB still holds. The major change is 
that the predicted contributions from BL Lacs and FS¬ 
RQs are unlikely to differ from zero. 

For the SFGs we have considered an alternative model 
in which the bias is set equal to that of a IO'^Mq dark 
matter halo. The bias of this object is larger than the 
reference value bspc = 1 at all redshifts. As a result the 
amplitude of the cross-correlation signal is expected to 
increase. However, the fractional increase is very small 
(< 10 %), and, due to the large error bars of observed 
cross-correlation data points, this change does not sig¬ 
nificantly affect our main conclusion that the IGRB is 
mainly produced by the SFGs. 

Our results seem to be consistent, within the uncertain¬ 
ties, with the outcome of different, independent analy¬ 
ses. Ajello et al. (2014) have been able to estimate the 
contribution of unresolved BL Lacs to the IGRB from 
their 7 -ray luminosity function measured from LAT 
data, and found that they do not account for more than 
10-15 % of the IGRB signal, consistent with our re¬ 
sults. In similar analyses focused on the FSRQs Dermer 
(2007), Inoue & Totani (2009), Inoue et al. (2010), and 
Ajello et al. (2012) have found that these objects pro¬ 
vide a similar contribution (~ 10 %) to the IGRB, again 


in agreement with our results, possibly increasing to 
~ 20 % when one accounts for objects with misaligned 
jets. 

As for SFGs, the estimate of Ackermann et al. (2012a) 
of a 4-23% contribution to the IGRB is consistent 
with our estimate of 20-95%, which, although favors 
a higher value, has a large uncertainty. A larger con¬ 
tribution from SFGs has also been recently suggested 
by Tamborraet al. (2014) which might be due to ac¬ 
counting for SFGs at z > 2 in the IR LF that have 
been recently observed by by Herschel PEP/HerMES 
(Gruppioniet al. 2013). It is also worth noticing that, 
within the 2 a error bar, our results are also consistent 
with those of Eields et al. (2010) that, based on the ex¬ 
trapolation of the 7 -ray production in the MW, find that 
SEG contribute to ~ 50 % of the IGRB (with rather 
large uncertainties). 

• In our analysis we have ignored MAGNs, even if they 
are likely to contribute the IGRB and its fluctuations, 
and restricted our modeling to SEG and blazars. The 
reason for this is that the expected cross-correlation sig¬ 
nal from these sources is robust to uncertainties in their 
bias parameters, whereas in the MAGN case model pre¬ 
dictions are much more sensitive to both their contri¬ 
bution to IGRB at high redshift and their (large) bias. 
Indeed, we find that, within the current uncertainties, 
their contribution to CCE and CAPS is degenerate with 
that of SEGs. One should keep this in mind when in¬ 
terpreting the results of our cross-correlation analysis. 
It might underestimate the expected cross-correlation 
signal at high redshift and, consequently, overestimate 
the SEG contribution whereas, in fact, part of the ob¬ 
served cross-correlation may be due to MAGNs. Pos¬ 
sible ways to isolate the contribution of MAGNs are 
more stringent observational or theoretical constraints 
on their bias and cross-correlation analyses with cata¬ 
logs of high redshift objects. 

The results of our work indicate possible directions for 
future research. Our analysis, which is mostly sensitive to 
sources at z < 2 suggests that the combined emission from 
SEGs, BL Lacs and FSRQs within this redshift does not com¬ 
pletely account for the whole diffuse IGRB signal. Extending 
our cross-correlation analysis to higher redshifts, using deeper 
catalogs of extragalactic sources can provide further informa¬ 
tion to clarify this scenario. 

While we observe a significant cross-correlation signal, the 
amplitude of the errors is still too large to efficiently dis¬ 
criminate among alternative IGRB models. We have learned 
that Fermi IGRB maps improve in both accuracy and preci¬ 
sion with time, not only because of the better photon statis¬ 
tics but also thanks to the revised Galactic Diffuse model, 
better characterization of the LAT PSF and to the identifica¬ 
tion and subtraction of an increasing number of point sources. 
We therefore expect that errors will be further reduced with 
the next Fermi data releases. Major improvements are also 
expected from multiwavelength catalogs, since the next few 
years will see the advent of next-generation galaxy redshift 
catalogs like eBOSS'\ DESI (Schlegelet al. 2011) and Eu¬ 
clid (Laureijs et al. 2011) extending over a large faction of 
the sky and containing tens of millions to billions of objects 
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with spectroscopic or photometric redshifts. With these fu¬ 
ture surveys, we not only expect to reduce the uncertainties in 
the cross-correlation analysis but also to be able to fully ex¬ 
ploit their tomographic potential which we have only started 
exploring in this work. 
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